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Stability Conditions in the Numerical 
Treatment of Parabolic Differential Equations 


1. Introduction. The numerical solution of a hyperbolic partial difference 
equation is in some circumstances subject to an “instability” which has been 
shown by Lewy and others! to have a simple significance with reference to 
the differential equation which the difference equation approximates. For 
example, the difference equation 


(1) (nina — 2bn,5 + bn 5—1)/ (At)? = (bags. — 2605 + bn-1,5)/ (Ax)? 
may be used to approximate the wave equation 
(2) ou = Coz. 


If initial values, ¢,,0 and ¢$,,1, are specified (corresponding to the specifica- 
tion of ¢(x,0) and ¢;(x,0)), equation (1) may be used to evaluate first ¢,, 2, 
then ¢,,3, etc. If a “component” of ¢,,;, i.e. a perturbation to ¢,,; which 
preserves the satisfaction of the difference equation and its boundary specifi- 
cations other than initial values, varies with m as (—1)* (i.e. with the short- 
est wavelength which can be represented with the lattice spacing Ax) and 
depends upon j as K‘, then from (1) 


(3) K? —-2K+1=-—40cK; o = (cAt/Ax)*. 


If ¢ > 1, equation (3) has one real root of magnitude greater than unity; 
hence this (shortest wavelength) component of ¢ grows exponentially with 
j. Thus if this component is minutely excited, say through rounding errors, its 
exponential growth may, for sufficiently large 7, destroy the resemblance of 
¢,,; to the solution of (2). The possibility of such an uncontrolled growth of 
error is termed ‘‘instability,”’ its absence “stability.”” For ¢ < 1 no such 
instability occurs. 

The significance of the abrupt change in behavior of (1) as cAt exceeds 
Ax may be seen as follows: Each ¢,, ; is computed by (1) with the use of pre- 
ceding values, ¢$n,;-1 and @$,,;-2, and the neighboring values, ¢,_:, ;-1 and 
$n+i,j-1- Thus ¢,,; is subject to influence only from within a “generating 
cone”’ described by j’ < j, |n — n’| < |j —j’|. It is ultimately determined 
by those initial values ¢,’,» for which |n’ — n| < j. In the difference equation 
the propagation of influences is thus limited to rates < Ax/ At. The differential 
equation (2) represents the propagation of influence along the characteristic 
lines, x + ct = constant; hence cannot be well represented by (1) unless 
Ax/At is at least as great as c. A more general study of hyperbolic difference 
equations shows the presence of instability whenever the mesh ratio Ax/At 
is too small to permit the propagation of influence along the characteristics 
of the corresponding differential equation. 

In the solution of parabolic difference equations stability conditions de- 
pending on the mesh intervals may also arise. They are considerably more 
burdensome than the above condition, usually requiring that At be small of 
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the order of (Ax).? By suitably modifying the difference equation, the limita- 
tion imposed by the stability condition can be removed. The character of the 
diffusion process described by a parabolic equation suggests that a stability 
condition which limits At/(Ax)? or even At/Ax is not fundamental. The 
effects of a localized disturbance are appreciable within a range which in- 
creases with elapsed time, ¢, about as #. Thus the region of appreciable in- 
fluence can asymptotically be included within an arbitrarily narrow gener- 
ating cone. 

In sections 2 and 3 the stability condition limiting the time interval in a 
conventional treatment of a simple parabolic differential equation is dis- 
played and its burdensome character illustrated. 

In section 4 an alternative simple parabolic difference equation approxi- 
mating the same differential equation is shown to be stable for all values of 
At/ (Ax). It is shown in section 5 that the condition for convergence to the 
solution of the corresponding parabolic differential equation is Ax — 0, 
At/Ax — 0. In sections 7 and 8 various generalizations of this method of 
calculation are examined, particularly with reference to their stability. In 
these generalizations the definition of ‘instability’ given above becomes 
inapplicable. The term will there be used in a looser sense to indicate a tend- 
ency toward the exaggeration of initially small ‘‘errors.’’ Accordingly the 
indications of stability given there are plausibility-proofs rather than 
rigorous arguments. 

2. Interval limitation for a one-dimensional diffusion equation. The 
diffusion equation 


(4) ¢: = P(x)bzz; p(x) > 0 
may be approximated by the difference equation 
(S) (Gnit1 — n,5)/At = Pa(bn—1,j —2bn,5 + bns,j)/Ax*?; Pa = p(ndx). 


If initial vaues ¢,,9 are specified, equation (5) permits the successive evalu- 

- ation of $n,1, ¢n,2, etc. An approximate stability condition for this process is 
readily derived by the assumption that p, changes slowly with m. An alter- 
nating component of @n, ;, dn,; ~ (— 1)", then appears in ¢,, ;,1 modified by 
the factor 1 — 4Atp,/Ax?*, 


On iti ~ (1 — 40pn)onj; o = At/Ax’*. 


An exponential increase of this alternating component may be expected if 
Pn > 1/20 over any appreciable region. The condition for stability is then 


(6) o < 1/2pm 


where pu approximates the maximum of the ,. More precisely, the pro- 
cedure described by (5) is unstable if an error in ¢,,; recurs in the same form 
in $»,j41 with increased magnitude, i.e. if an error component 


5¢nj = Kn; |K| >1 


satisfies (5) and \,, satisfies homogeneous boundary conditions corresponding 
to those imposed on ¢,,;. For definiteness it may be assumed that @o,; and 
ow, are specified. Then an “‘error-eigenfunction,” \,, must satisfy 


(7) = ADrn = (K = 1)A,/o = PalAn—1 - 2An + An+1)3 Xo = Aw = (. 
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Multiplying by A,//, and summing yields 


N N 
~— 4p pa An’/Pn ie Drn(Angs ne An) mel An(An — An-1) ] 


n=0 n=0 
N-1 N-1 
” x An(An41 _ An) = } Anti(Angi = An) 
n=0 n=0 


N-1 
= x (Anti nal An)? < 0; 
n=0 
thus the positiveness of an eigenvalue, ~, follows from the assumed positive- 
ness of p,. This requires that K < 1; hence instability can arise only if for 
some error-eigenfunction K <— 1. The stability condition (6) holds 
strictly if ps denotes the greatest eigenvalue, p, of (7). 

The stability condition (6) imposes a burdensome limitation on the 
numerical treatment of parabolic differential equations by this method since 
a moderately fine spatial division may require an immoderately fine temporal 
division. This will particularly be the case if a very fine structural repre- 
sentation is required (large N) or if a great disparity exists between the 
largest and smallest values of p,. If the range in ¢ for which a solution is re- 
quired is of the order of the time required for effective diffusion to the dis- 
tance NAx, the number of time intervals required will be of the order of 


tmax/At ~ N*pu/ Po 


where $o approximates the smallest of the pn. 
3. Application to the radial diffusion equation. To illustrate the above 
general conditions we consider the simple diffusion equation 


(8) b: = der + ¢,/7; $(r,0) = &(r) 
o(1,t) =a 
o(R,t) = b. 
Writing r = e* puts (8) in the form 
(9) ot = Chez; $(x,0) - (x) 
¢(0,t) =a 
o(InR,t) = b. 


This is of the form of equation (4) with p(x) = e~**. The error-eigen- 
function equation (7) is now 


(10) 4pnn _ (1 as K)n,/¢ > Catal ( ee + 2n. + Nn—1); no = nv = 0, 


where 9, = (— 1)"Aq. 
The greatest eigenvalue, pu, is readily approximated by a variation 
principle : 


For each solution of (10) the expression 
N-1 N 
x=2 a na(ta + natt)/D oe" 
n= 0 


takes on a stationary value, namely 4). 
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The maximum value permitted to x is then 

(11) Xu = 4pm = (1 — Kmin)/o. 

The condition for stability, Kmin > — 1, then requires 
o£ 2/xy = Cerit- 


A lower bound to x,,, hence an upper bound to geri (for N = ©) has been 
obtained by setting 7, = me~*" and maximizing x with respect to 6. The 
resulting approximation to ocrit is displayed (vs. Ax) in Fig. 1. To estimate 
the error in this approximation numerical solutions of the difference equa- 
tions for two values of Ax were carried out with excessive o’s. The rates of 
growth of the dominant instabilities determine x, (by eq. 11), hence 
Gerit- These two ‘“‘experimental’’ values are also shown in Fig. 1. 
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’ To illustrate the uncomfortable limitation imposed by the stability 
condition we assume (arbitrarily) that the calculation is to be carried to 
sufficiently great ¢ that ¢(7,t) at r = 100 is substantially affected by the 
inner boundary value, a. This requires tmax ~ (100)*. To represent ade- 
quately the structure of ¢ near r = 1 requires, say, Ax = 0.1. Then the 
stability condition requires 


At < oerit(Ax)? = 0.0085. 


The number of ‘‘cycles’’ (j-values) required for the calculation is then 
tmax/At ~ 10°. A moderately fast electronic digital computer could complete 
the calculation in a few hundred hours! 

4. A generally stable difference equation. The difference equation (5) 
determines a set of numbers, ¢,;, which approximate the solution of the 
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differential equation (4) at the junction points of the rectangular lattice 
x = nAx,t = jdt. The alternative difference equation here considered makes 
use of a diagonal lattice, obtained by omitting those junction points for 
which n + j is, say, odd. The term ¢; is now represented by the difference of 
two ¢,;-values of the same m and j’s differing by 2; ie. o:(mAx, jAt) ~ 
(bn, 41 — On,j-1)/2A4t; m + 7 odd. In the representation of ¢22 by a second 
difference with respect to m the end terms, ¢,_1,; and ¢.41 ;, may be used as 
previously but the term, ¢,,;, corresponds to a point omitted from the lat- 
tice. It is therefore replaced by the mean of the two terms of neighboring 
j-values. 


dez(nAx, jAt) ~ (bn-1,5 — On j—-1 — Onig1 + On41,5)/(Ax)?; m+ 7 odd. 


The resulting ‘difference approximation to (4) is 


(12) Onset — On ja = 2hno(Gn-1.5 — On j—1 — On sti + ners), 
n+jodd, p,>0, o = At/Ax? 
which may be written as 


(13) Gn, i41 = On j—1 + On(On—-15 — 2Gn, 5-1 + Ont13)» 


where ; 
On = 2pnc/(1 + 2pac); 0 < a, < 1: 


If the differential equation has terms in ¢, and ¢ these may be represented 
by (¢n41,5 — On—i,s)/2Ax and ($n, 541 + On, j-1)/2 respectively. 

It is to be noted that the initial conditions required to permit calculation 
with (13) require the specification of ¢-values for two initial times (j-values), 
or initial boundaries, x(#), when only one would be required by the conven- 
tional difference equation. This requirement suggests that (13) is of hyper- 
bolic rather than parabolic character (cf. section 5). The second set of initial 
values may be computed from the first set by use of the conventional 
equation with sufficiently small At. 

In each cycle of the calculation ¢,,; may be evaluated for all n-values 
associated with one j-value, even m occurring in one cycle, odd m in the next. 
Alternatively the ¢,,; may be evaluated in each cycle for all m along a 
diagonal line » + 7 = constant (or m — j = constant). These two orderings 
of the calculation, called the “leap-frog’’ and “‘pyramid’’ methods respec- 
ively, are shown in Fig. 2. 

In the pyramid method the cycle number may be used to index the 
¢-values in place of j. If the cycle number, m, is written as a super-script, 
the difference equation (13) takes the form 


(14) dnt! = oe + an(Gri — 2m + ory). 


This form of the difference equation displays clearly the similarity to the 
“extrapolated Liebmann”’ form of the relaxation method.? 

The qualitative behavior of the solution of (13) may be investigated by 
examining the propagation of a spatially sinusoidal component of ¢, treating 
Pn, hence also a», as a constant. If 


Ags = Kiet» 
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satisfies (13) for constant a, then 


K = K7 4+ y(e-” — 2K— + e*). 
Hence 


(15) K? — 2Kacos a+ (2a — 1) = 0. 


For each y, two values of K satisfy (15). Thus a sinusoidal component of 
¢ has two modes of propagation with increasing t. The propagation factors, 
K, are shown in Fig. 3. For y near 0 or near x the propagation factors are 


cycle # 
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real, for intermediate y they are complex. In no case is the magnitude of K 
greater than unity. The difference equation (13) is thus seen to be stable for 
any constant positive value of px. 

To show the stability of (13) for a varying a,, we assume the presence 
of an error-eigenfunction of the form 


oni = Kin; Xo = Aw = O. 
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Substitution in (13) yields 


Ae-1 + Aa(i _ K? = 2an)/Kay + Anti = 0 
or 


(16) w?(2An - An-1 7 An+1) —_- 4w(az* _ 1)A, 
+ (2rn + Ani + Anti) = 0, 
where 
(17) K = (w — 1)/@ + 1). 

Multiplying (16) by A,* and summing over m = 1,2,---, N — 1 pro- 
duces, after slight rearrangement of terms, the equation 


N N-1 N 
(18) wD [An — Anal? — 40 DL (ae* — 1)/An|*? + D [An + Anil? = 0. 
1 1 1 
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The positiveness of the three summations ensures that the real part of » is 
positive, hence by (17) 


(19) |K| <1. 


If the fixed end conditions, A» = A, = 0, are replaced by the homogene- 
ous boundary conditions, 


Ai — Ao = Ado me 
aes issih csieies A and B real and positive, 


equation (18) is modified by the inclusion of additional positive terms in the 
first and third summations. Then again (19) follows from the positiveness of 
(a,;' — 1); hence the difference equation (13) is stable for any positive co. 
5. Dependence of the solution on the lattice spacings. To illustrate the 
dependence of the difference equation (13) and its solution on the lattice 
intervals, Ax and Af, we examine in detail the special case of constant p(x). 
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This constant may be given the value unity without loss of generality. The 
difference equation is then 


(20) dn sti — Gnj—1 = 20 (bn-1.5 — bn j—-1 — bn iti + baits), ¢ = At/Ax’, 
which approximates the homogeneous, one-dimensional diffusion equation 
(21) dt = Grz. 


If the set of four ¢-values connected by (20) are assumed to be imbedded in 
a function ¢(x,t) which permits a power series expansion about the ‘‘center 
point” (m,j), they may be described by the expansions 


1 i 1 
Pnitj = o+ Ax: + 2 (Ax)*z2 = 6 (Ax) *z22 + 24 (Ax)*dzz2z +o, 


1 1 1 
on, jz1 = & + Albi + 3 (At)*oee ra (At)*peee + 24 (At)*beeee + °° 


Substitution into (20) then yields 
1 1 
(22) ot — ozz = 12 (Ax)*ozzz2 id (At)*be/ (Ax)? aie 6 (Ab) bece 


1 
— 2 (At)*preee/ (Ax)? +-->. 


The use of (20) as an approximation to the diffusion equation thus con- 
sists in the neglect of the right member of (22). The use of small values of 
Ax and At is not alone sufficient to justify this neglect. To make the second 
term of the right member of (22) negligible requires also that the ratio, 
At/Ax be small. If Ax and At approach zero with constant ratio, c = Ax/At, 
the difference equation (20) does not approximate the parabolic equation 
(21), but rather the hyperbolic equation 


(23) u/c + bt = daz. 


It may be noted that the mesh ratio, At/Ax, is just the maximum permitted 
by the stability condition for the treatment of this hyperbolic equation by 
the usual difference method. The stability of the difference equation (20) 
may be regarded as arising from the introduction of this hyperbolic term in 
the approximating differential equation. 

The above considerations indicate that the requirement for accuracy of 
treatment of the parabolic equation approximated by (20) is the smallness 
of both Ax and At/Ax. The ratio, At/Ax is not, however, required to be small 
in the order of Ax, as would be needed for a fixed value of o. Since (21) 
implies 

bee = dzezz 


the dominant terms in the right member of (22) may be written 
1 
bt — bes = 75 (AX) *basee(1 — 120%) +---. 


Thus for fixed Ax the optimum value for At, from considerations of accuracy 
alone, lies near 12-4Ax*. 
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The accuracy of the difference equation representation of (21) may be 
examined by reference to the Fourier expansion used in establishing its 
stability. A component of @ having the sinusoidal dependence, e**’, is 
exponentially attenuated with increasing j by the factor K per cycle. The 
attenuation factor, K, is determined by (15) and is double valued (again 
suggesting the hyperbolic character of the difference equation). As shown in 
Fig. 4, K,, the greater value of K, closely approximates the correct attenu- 
ation factor, K. = e~”** (i.e. the factor by which a sinusoidal component 
e*77/4z is attenuated in time At by the differential equation (21)) for small +. 
With increasing y this greater K becomes progressively more unrealistic, 


te) 





005 4 


1 4 


“s+ 














Fic. 4. 


the more rapidly the greater is ¢. For 
0 < cosy < (1 — 1/40*)3 


the,two K’s become complex, of common magnitude [(2¢ — 1)/(2¢ + 1) }!. 
For negative values of cos y the attenuation factors are the negatives of 
those for the corresponding positive cos . Since only even values of m + j 
are considered these negative values of cos y are superrogatory. 

The rapid increase of the error in K, with increasing y would seem to 
indicate that the difference equation (20) is of unacceptable accuracy except 
for small o. After a large number of cycles of calculation (large 7), however, 
the two values, K, and K,, appear in the solutions of (20) and (21) raised to 
the power j. For sufficiently large j, both K,/ and K,/ are essentially zero 
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except for very small y. An approximate measure of error is thus the greatest 
value of the discrepancy, 


D; = Ki _- K,/ = e7 is? _— K,!. 
The solution of (15) in power series in y* yields 
K, = K.exp[— y*(o* — o/12) +---]. 


For 7 > 4(¢ — 1/12c) the greatest discrepancy occurs at y? = 2/jo and 
has the value 
D; = 4e?(¢ — 1/12¢)/j. 


The difference equation (20) thus permits increasingly accurate representa- 
tion of the solution of the differential equation (21) as 7 becomes large in 
comparison with 4c. 

The above estimate of accuracy takes into account only the errors in the 
greater attenuation factor K,. An additional error may arise from the fail- 
ure of the initial specification of ¢,, ;-values to associate the correct amplitude 
with the components which decay with attenuation factors, K,. Since 
initial ¢-values must be specified for two successive j-lines it is not sufficient 
that the ¢’s closely represent the initial values of the corresponding para- 
bolic differential equation. It is also necessary that the initial time deriva- 
tive be well represented. The latter requirement is the more stringent the 
larger is o, since K, and K_ are then closer together, hence more difficult to 
discriminate in the initial-value specification. The value ¢ = } (hence 
a = 3) is especially favorable in this regard, since for it K_ is identically 
zero (cf. eq. (15)). The special property of a = 3 is also to be seen in (13). 
For a = 3 the terms in ¢,, ;-1 drop out, hence specification of initial ¢-values 
for one j-line suffices. 

Apart from the special case of constant p,, no choice of ¢ makes a, uni- 
formly equal to 3. It nevertheless seems expedient to begin the calculation 
with a o which makes 3 a typical value for a,. After a sufficient number of 
cycles ¢ can conveniently (i.e. without requiring interpolation) be increased 
by an odd integral factor. 

6. Examples of the use of the generally stable difference equation. The 
first example chosen to illustrate the use of the difference equation (13) is 
the equation of diffusion to the boundary of a one-dimensional semi-infinite 
medium. Its differential equation and boundary conditions are 


¢: = dz fort>0, O<x< ~, 
¢(x,0) = 1, $(0,t) = 0. 
It has the simple solution 
(x,t) = erf (x/2t4). 


Numerical solutions were obtained by the pyramid method for two 
values of o, $ and 5. The boundary values are 


(24) oo" = 0, ¢,° = 1 forn>O. 
Subsequent ¢-values are determined by the relation 
(25) ont? = on + an(Grpi — 2dr + dr-1), 


a, = tor 10/il aso = $orS. 
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test For each cycle (m-value) the ¢7*' are computed for successively smaller n, 
beginning with an n sufficiently large that ont does not differ appreciable 
from unity. (Here 2m = n + j.) 

The resulting ¢-values are displayed in Fig. 5 and 6, plotted for several 
values of t = At(2m — n) vs. the similarity variable, 


p = x/2th = n/2[o(2m — n) }. 


and The approach of ¢,” to ¢(u) with increasing 7 = 2m —n, hence a 
fortiori with increasing t = ojAx*, is noticeably more rapid for the smaller c. 
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Fic. 5. Results of One Dimensional Diffusion Problem 
A second example used is the radial diffusion equation, 
wo (26) ot = der + ¢,/7; r>1. 
Setting r = e* permits writing (26) in the form 
(27) o: = bss; x > 0. 


The boundary values 
$(0,t) = 0, o(x,0) = 1 
were chosen. 
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Again the pyramid ordering was used; hence the difference equation is 
(25) with 
Qn = 2ce?4*/(1 + 2ce-2"4*), 


The boundary conditions of the difference equation are given again by (24). 
Three solutions were obtained, with intervals as follows: 


Case Ax At o = At/Ax* 
I 05 .00125 i 
II 05 .0125 5 
III 10 05 5 


The solution in Case I shows good accuracy even in the first few cycles. 
A comparison with an analytically derived solution for ¢ = 1/100 = 8At is 
shown in Fig. 7. The four points plotted are obtained in the fifth, sixth, 
seventh, and eighth cycles respectively. 

A somewhat greater number of cycles is required for comparable ac- 
curacy in Case II. 
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Its solution is shown for ¢ = 5A¢ and ¢ = 17At in Fig. 8. For the latter 
time the plotted points are obtained in the ninth and subsequent cycles, and 
are of about the same accuracy as the points in Fig. 7. 

The solution in Case III is shown in Fig. 9 for ¢ = 5A#, 15A#, and 40Az. 
The values for 40A¢ = 2.0 (derived from the twenty-first and subsequent 
cycles) are not distinguishable in this graph from the previous cases nor 
from the analytic solution. 

7. Extension to other parabolic equations. The demonstration given 
above of the stability of the leap-frog (or pyramid) method depends upon 
an eigenfunction expansion of the errors in the dependent variable. (So, in 
fact, does the definition of stability there used.) To justify the use of these 
methods for broader classes of parabolic equations other indications of 
stability must be sought. In section 8 several types of linear parabolic 
difference equations are shown to display properties indicative of stability. 
No uniform definition and proof of stability has been produced, however. 
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A similar computational technique can be applied to a nonlinear equa- 
tion, 


(28) o: = f(¢, gz, $zz, X, t). 

In suitable cases ¢(x,t) may be sufficiently well represented by a set of num- 

bers, ¢,,;, satisfying the corresponding difference equation, 

(29) (bn, i41 — bn,j-1)/2At = fL3 (bn, i41 + bn5-1), 3(bn41,5 — bn-1,5)/Ax, 
(On41,5 — Onj-1 — Oni41 + On—1,5)/(Ax)?, ndx, jdt]. 


The stability of this computational procedure depends upon the rates of 
growth of small deviations from this ‘‘correct” solution. If equation (29) 
permits a power series development in these deviations then they will, when 
sufficiently small, be governed by a linear parabolic difference equation of 
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the same form. The above stability arguments for linear equations thus 
indicate a wide range of usefulness for this procedure as applied to nonlinear 
equations. 

To illustrate this procedure we consider the equation describing the one- 
dimensional isothermal flow of a perfect gas in a porous medium. 


(30) ¢: = (@)ez for t>0, x>0 
o(x,0) = 1; $(0,t) = 0. 
A convenient corresponding difference equation is 
Ad = On, 541 — On,j-1 = 2o[P naj — 26? 0 j-1 + Pntt js — Ab(bn-15 
+ Pant L i]. 
The use of $n-1,; + Gn4i,; rather than oq, ;-1 + $2,441 as the coefficient of 
A¢ in the right member is consistent with the order of approximation to (30) 


4 4 4 i 4 4 
t T T ? * t 
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with successively increased 
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used and simplifies the explicit expression for @», j+:. This is 


(31) dnjna = bn, j—-1 + 2o[G n-1.5 — 26? a j-1 + nts] 
/(1 + o(@n-1,5 + n41,5)]. 


Solutions to (31) were obtained with o initially given the value 1, later in- 
creased by successive factors of five. (An attempt to start a solution with 
a = 5 led to uncontrolled oscillation.) When o was not increased too soon, 
qualitatively correct solutions were obtained. The solution of (30) is of the 
form g(x/t#). Accordingly, to provide a measure of accuracy of the approxi- 
mation of ¢,,; to the solution of (30), the ¢,,; values for each cycle were 
used to determine a graphically interpolated value of § = /j* corresponding 
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to the arbitrarily selected value, ¢,,; = 0.6. The discrepancy between ¢ and 
its asymptote, approximately 0.586, indicates the order of inaccuracy of the 
solution. Fig. 10 shows values of ~ for a family of solutions of (31) differing 
in the number of cycles carried out at each value of «. The abscissa, dis- 
played logarithmically, shows ¢ in units of the initial At rather than the 
number of cycles of calculation performed. The error (of approximation to 
the solution of (30)) is seen to increase and later subside slowly after each 
increase in o. The peak error decreases rapidly with deferral of the increase 
in o. To keep the error of the order of one percent requires a few hundred 
cycles at each o-value, somewhat more in the later stages. 

8. Stability for the general parabolic linear equation. The differential 
equation, 


(32) | bt = P(x) bez + a(x,t)be + 1 (x,t) p 
may be approximated by the difference equation 


(33) (bn jt1 — On j—1)/2At = Paz (Ong1.5 — Onjt1 — Onj—-1 + Gn—15)/ Ax? 
+ Qnj(bno1,5 — Gn—1,4)/2AX + 125 (On, i41 + On5-1)/2; m+ 7 odd, 


where the notation and lattice structure is as described above. The solution 
of (33) may be expected to provide a reasonable approximation to that of 
(32) only if throughout the region considered 


(34) p20 
(35) |gAt/Ax|< 1, and 
(36) \rAt| K 1. 


It is assumed here that initial rather than final values of ¢ are specified ; 
hence (34) is necessary for the uniqueness of the solution. Condition (35) is 
clearly required by considerations of causality similar to those described 
above. (i.e., the lattice slope, |Ax/Az|, must clearly be at least as great as 
the ‘“‘stream velocity,” |g|.) It ensures that the cone of determination of 
each ¢,; includes the characteristic curve of (32), satisfying 


dx/dt = — q(x,t) 


and passing through (nAx, jAt). In the event p = gq = 0, (32) becomes an 
ordinary differential equation (involving x as a parameter) which is reason- 
ably approximated by (33) only if (36) is satisfied. 

We first consider the case g = r = 0. Equation (33) then reduces to 
(37) (Gn, i+1 — Oni) = (bn, 5-1 — Gas) (1 — 2Paso)/(1 + 2pn;0) 


where 


(38) Onj = 3(bn—-1,5 + Ony1,5) and o = At/Ax? 
hence 
(39) lon, s+1 — Ons| < [bn s-1 — Pail. 


Each stage of the computation is characterized by a “front” in the n,j 
lattice bounding the region for which ¢-values have been computed. The 
front is specified by a single-valued function, j(m), (7 + even), satisfying 
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j(n +1) =j(n) + 1. For j < j(n), oa; has been computed, otherwise not. In 
each step of the caiculation a “low”’ point of the front, j(m) = j(m + 1) — 1, 
is replaced by j(m) + 2. Thus a “‘valley” in the front is replaced by a “‘hill.” 
We define a function of the front consisting of the sum of squares of the first 
difference of @, taken along the front: 


(40) F = D> (bnin) — On—1, (2-0). 


In a step of the calculation two terms of F are altered as a result of the re- 
placement of ¢a, jn) Dy a, i¢m42. With the use of (39) it can be shown that 
this replacement diminishes F unless p,; = 0 or n,j41 = $,j, in which 
event F is unchanged. Thus as the calculation proceeds F decreases mon- 
otonically except by reason of influences introduced at the boundaries. This 
non-increasing character of the (positive) function, F, may be regarded as 
indicating the stability of (33). (Here, and below, the term “stability” is 
used loosely to indicate the absence of any tendency toward disastrous in- 
crease of irregularities in ¢,,; initiated by rounding errors and the like. It is 
hoped that, pending development of rigorous criteria, these imprecise indi- 
cations of “stability in the mean” may provide useful hints to the working 
computer.) 
A similar indication of stability can be displayed for the case, 


(41) |gAt/Ax| £1; |gqAt/Ax| < 2po; r= 0. 
This permits writing (33) in the form: 

(42) (1 + 2po)y4 = (1 — 2po)y_ + gAt/Ax 
where 


(43) V+ 7 (bn, 541 oY $n, i)/(On+1.5 a On-1,4)3 y- a (On, 5-1 a $n, i)/ 
(Gn+1,4 — On—1,5)- 
Thus either 


(44) W+] < WW] or ls] < 4. 


This permits the inference that |@n, j;-1 — ¢n—-1,5| + |@n41,3 — Onj-1| is (in 
the first event) decreased or (in the second event) at least not increased as a 
result of the replacement of $n, ;-1 by @n,;41. We can thus define a positive 
front-function 


(45) G= >> In, in) ous On—1, (n—0| 


which is non-increasing (except by reason of disturbance at the boundaries). 
The operation of the difference equation thus tends to bring the ¢-distribu- 
tion on the front to minimum total variation. This property gives a weaker 
indication of stability than that described above since G, unlike F, is not 
typically minimized by a unique ¢-value. It nevertheless seems to justify 
describing the difference equation as stable. 

The second restriction of (41) limits the usefulness of this argument if 
2po becomes small. An argument which is then applicable may be made as 
follows: Define 


(46) Mt = Onji-1 — On-153 Mar = Onit5 — Onj-1 
Prt = onj41 — On-1,53 Pr = On41,5 — Onj41- 
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Then equation (33), with r = 0, may be shown to imply 
(47) (1 + gAt/Ax)Pp? + (1 — gAt/Ax)Py? = (1 + gAt/Ax)M 7? 
+ (1 — gAt/Ax)Mz? — 4po(1 + 2pc)2{ Mz + Mr 
+ (Mi — Mpr)qAt/Ax}?. 
A front-function which displays a decreasing tendency may thus be written 
as 


(48) H => (1 + gAt/Ax) (bn, jn) — On—1, (nv)? 


the + asj(m — 1) = j(m) +1. 


In (47) g is properly written q,;. In (48) it becomes ambiguous, being either 
Qn—1, i(n) OF Yn, j(n)41- Thus we can only assert that H is non-decreasing except 
by reason of variations in g with the advance of the front or disturbance in- 
troduced at the boundaries. 

In equation (33) with nonvanishing 7, the precise meaning of ‘‘stability”’ 
is not evident and no demonstration of properties approximating this con- 
cept are known to us. A preliminary qualitative examination failed to dis- 
close any indication that variations in ¢ can grow to an alarming extent. 
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On Over and Under Relaxation in the Theory 
of the Cyclic Single Step Iteration 


In order to speed up the sometimes very tedious computations in solving 
equations by the single step method, the device of the so-called ‘‘incomplete 
relaxation” (under or over relaxation) has been often used, although appar- 
ently no systematic discussion of this device has been as yet tried." 

It may seem, however, that in the case of the “relaxation procedure’’ the 
speeding up can be achieved in this way only in special cases, at least in the 
case of a symmetric positive definite matrix. Indeed, if the progress of the 
computation is measured by the decrease of a corresponding quadratic form 
A (g%) depending on the k-th approximating vector ¢;%, we have the formula 


(1) A (Se) — A (Sess) = ge(2 — ge) |rP.)*/an,w, 
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where JN, is the index of the variable the value of which is improved at the 
k-th step, 7, the corresponding “‘residual”’, g a coefficient characterizing the 
degree of the incomplete relaxation at this step. (We have usually 0 < q& 
< 2; for 0 < qe <1, we have under relaxation and for 1 < q@ <2 over 
relaxation, while for gq = 1 we have the usual complete relaxation.)? 
However, the above formula shows that the decrease of A({x) is maximum 
if qx is taken as 1. Thus far the improvement in the case of gq, * 1 appears 
to be only possible if for such a value of g: one of the residuals at the next steps 
will come out particularly large. 

In what follows we discuss completely the case of a real 2 K 2 matrix, 
symmetric or not. It then turns out that, if the cyclic single step iteration 
converges, this convergence can indeed in all cases be essentially improved 
in using the incomplete relaxation with a convenient coefficient g, the same 
at each step. It turns out that even in many cases of unsymmetric matrices 
where the usual cyclic single step iteration diverges, it is made convergent 
in using a convenient incomplete relaxation. 

The usual cyclic single step iteration is applied to a system of equations 


(2) L Gy X = Vu (u = 1,---, m) 
v=] 


with a non-singular matrix A = (a,,). We can assume without loss of gen- 
erality that y, = 0 (u = 1,---, 2), since we can always obtain this by chang- 


ing the origin. Then we get from one approximation vector § = (x1,---, Xn) 
the next one &’ = (x;’,---, Xn’) by the formula 

n p—1 
(3) Onn Xp = Cuyp Xp — LD Ope X» — L Spr Xo’ (wu =1,---, n). 

v= v=1 


This is replaced for the incomplete relaxation with a coefficient g by 


n p-l 
(4) Cup Xp = Guy Xp — J Dd Gur Xr — J DL Aye Xr’ (u = 1,---, 2) 
=p v=! 


and this can be written in the form 


ul 
(5) yp Xp + Lo Oy Xp’ = Ay (1 — Q)%, — A Dyy Xy. 
v=1 vamp tl 


We now decompose A into the sum 
(6) A=L+D+R 


where D is the diagonal matrix containing the main diagonal of A, while in 
L all elements on the diagonal and to the right of it and in R all elements on 
the diagonal and to the left of it vanish. We make here the assumption, 
which is usually made in the theory of the single step iteration, that none of 
the diagonal elements of A vanishes; then we can write (5) in the form 


(D + gL)?’ = (1 —g@D — QR) 


and in solving this, since the triangular matrix D + gL is non-singular, we 
obtain 
t = (D+ ql)“ ((1 — gD — @RE. 
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If now we put 
(7) Q. = (D + gL) ((1 — g)D — gR) 
we see that the k-th iterated vector & is obtained from the starting vector 
& by the formula 
& = Q,*to (k = 1,2,---). 

In order that our iteration be convergent for any starting vector £o, it is 
necessary and sufficient that the maximum modulus A, of the characteristics 
roots of the matrix Q, be less than 1. The speed of the convergence is then 
measured by A,. The smaller A,, the faster is the convergence. 

The characteristic equation of (7), 

|AE — (D + gL) ((1 — g)D — gR)| = 0, 
becomes, if the matrix of the left hand expression is multiplied by D + gL, 
(8) |A(D + gL) — ((1 — g)D — gR)| = 0. 

We discuss this equation only in the case of a real matrix A with n = 2. 
Here the number 
(9) & = (@42021)/ (@11422) 
is the characteristic constant of the problem. Our results are contained in the 
following four theorems. 


THEOREM 1. Jf u > 1 then for all q from (0,2) we have A, > 1 and the 
process is divergent. 


THEOREM 2. Suppose that u <1 and put 
(10) go = 2/(1 + (1 — u))). 


Then with monotonically increasing q, Ag is monotonically decreasing for 
qg < qo and monotonically increasing for q > qo, so that the optimal value of 


A, is 
(11) Aopt = Ag, = |u|/(1 + (1 — u)!)? (u < 1). 
THEOREM 3. Suppose that u < 0. Then a necessary and sufficient condition 
for the convergence of our procedure is that 
(12) 0<g<gqi = 2(1 + [uli 


If 0 < u < 1, we have convergence for all q with 0 < q < 2. 

THEOREM 4. Suppose that |u| <1. Then a necessary and sufficient condition 
for Aq < Ai ts that q be contained in the interior of the interval between 
land1i+u. 


To prove these theorems we start from the corresponding case of (8) 
(A+ . one 1)au qai2 

Agao1 (A + — 1)a22 
Without loss of generality we can assume that a1; > 0, do. > 0. In 


dividing the first row and the first column by a;;! and the second row and 
second column by d2:! we can reduce A to the form 


a=(° > u = af. 
a 1 














(1 


ar 
tr. 


2 iQ ¢) © -os 


~~ me 0 oe 


tor 


it is 
tics 
hen 


qL, 


the 


the 


for 


e of 


lion 


and 








CYCLIC SINGLE STEP ITERATION 155 


Our equation for \ becomes now 
A+q-1 gB 
Aga A+q-1 
= + A(2g — 2 — gu) + (¢ — 1)? = 0. 


(13) N,(A) _ 


For g = 1 we obtain 
(14) Ai = |14| 


and we see that the usual cyclic single step iteration converges, for an arbi- 
trary starting vector, only if |u| < 1. 

In the following discussion we assume u ~ 0. The discriminant A of the 
polynomial N,(A) is 


1 
A =4¢u (ug? — 4q + 4). 


_ This vanishes for g = go = 2/(1 + (1 — u)}), where obviously 


(15) Sept (0 <u <1) 


0<q <1 (u <0), 


while the other root 2/(1 — (1 — u)#) exceeds 2 if 0 < u < 1 and is negative 
if u < 0. We see therefore that - 


| A>0O (u > 1) 
Sgn A = Sgn u(go — g) (u <1). 

We consider now two cases according as A < 0 or A > O. By (16), Ais 
negative only if « < 1 and either 0 < u < 1 and q@ <q <2 o0ru <0 and 


0 <q <q. In both cases we have obviously from (13), A, = |¢ —1| 
(A < 0) and in particular 


oe ab (u>0, 2>q> 4q) 
A,=1-q (u<0, 0<@< q). 


In virtue of (16), the hypothesis u > 1 of theorem 1 is never realized for 
a negative A. Theorem 2 follows immediately from formulas (17). Theorem 
3 follows from the fact that by (17) A, < 1, while for u < 0, in any case 
go < gq: and g cannot exceed go by (17). It follows finally from (17) that for 
u > 0 the condition 


(16) 


(17) 


is equivalent to g < 1 + u while g remains > go > 1. On the other hand, if 
we have u < 0, the condition (18) is equivalent to g > 1 + u, while g re- 
mains <qo < 1. Therefore theorem 4 and all our assertions are true in case 
A<0. 

Now let A > 0. If \ is the root of N,(A) with |A] = A,, we have 





“6 
(19) h=R+ <A, R=S(¢- 245 ). 


where « = Sgn R. We have obviously 
(20) Ag = ea. 
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The monotonicity of \ with respect to g depends on the sign of 


dd _ , qua — »— (q — 1) 
dq ~2k—qu—2+2q° 





By (19), the denominator is equal to 
2(. — R) = 2A} 


and has the sign of «. If we denote the numerator by 6, we have 


(21) §=A(qu —1)+1-q 
and therefore, by (20), 
(22) Sg z = Sgn 6. 


We deal first with the case u > 1. Here we have A > 0, the roots of 
(13) are real and, since 


N,(1) =1—qu—-2+2¢+¢@-—2¢+1=¢(1 —u) <0, 


one root A exceeds 1, and so A, > 1. We see that in this case we have diver- 
gence for any value of g > 0 and theorem 1 is proved. 

Since for u = 1, A becomes singular, from now on we can make the 
assumptions 


(23) u<i1, gqoisreal, A>O. 
We have then in particular from (16) 


either u>0O, g<qo, go>1 


(24) or 4u<0, g>qQ, <1. 


We prove now the following lemma, the proof of which is the main 
difficulty of the paper: 

Lemma. Under the assumptions (23) ué is always negative. 

Denote by 6; and 62 the two values of 6 corresponding to the two roots of 
N,(A). We have from (21) and (13) after some simplifications 


(25) 51 + 62 = qu(qu — 3g + 2), 
(26) 5:62 = (q — 1)uqg?(1 — x). 

The expression on the right in (25) vanishes for g = g2, where 
(27) gz = 4/(3 + (9 — 8u)#) <1 (u <1), 


while the other root exceeds 2 for u > 0 and is negative for u < 0. Therefore 
we have 


(28) Sgn (6; + 62) = Sgn u(g2 — Qq). 


We will now consider separately the cases u > 0 and u < 0, and prove 
in the first case 6 < 0 and in the second 6 > 0. 

In the case u > 0 we have 1 > u > 0, and, since A > 0, g < go. From 
(24) and (27) it follows that 


(29) 2>q>1> qQ. 
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If gq exceeds 1, it follows from (29), (28), and (26) that 
6+6.<0, 66:>0, 5<90, &<0 


and therefore 6 < 0. 

If g < 1, it follows from (26) that 6,5. < 0. On the other hand the ex- 
pression of R in (19), for u > 0, g < 1, becomes positive. We have therefore 
in this case « = 1, A is the greater root of N,(A) and since gu — 1 in (21) 
becomes negative, we have 


§= Min (61, 52). 


But then it follows from 6:6. < 0 that 6 is negative in this case also. 

We consider now the case u < 0. Since A is assumed positive, we have 
here by (16), g > qo, while on the other hand from (24) and (27) it follows 
that 


Gz <q <1. 
Since therefore in any case g > go, we have from (28) 
(30) 8, + 52> 0. 
On the other hand we conclude from (26) 
(31) Sgn 6:52 = Sgn(i — gq). 


If therefore g < 1, we have 5:5. > 0 and from (30) 
5,>0, &>0, 6>0. 


Suppose now g > 1; since in this case, by (19), R is negative, we have 
e =— 1, Ag =— \ and from (21) it follows now, since the coefficient of \ 
in (21) is negative, that 
$= Max(éi, 52). 


But now we see from (30) that 6 is again positive. Hence our lemma is 
proved. 

We see now, from (22), that in the case of positive A, A, monotonically 
decreases for g < go and monotonically increases for g > go, as we already 
deduced from (17) in the case of negative A. 

The minimum of A, is obtained for g = go. We have, in applying for 
instance (17) and in using (10), 


(32) Ag, = |(1 — (1 — u)4)/(1 + (1 — u)4)| 
= |ul(1 + (1 — u)4)* (u < 1), 
and theorem 2 is completely proved. 

The expression. (32) is less than |u| unless wu = 1. For |u| < 1, we there- 
fore always have an improvement in using incomplete relaxation with 
@ = %, which is particularly pronounced for small |u|, since we have 
(33) Ag ~ |u|/4 (u 0). 


As to the values of g for which we have convergence at all, we have, since 
Ao = 1, convergence when q is in the interval (0, go). On the other haud, 
the roots of (13) become, for g = 2, 


(34). 2u—1+2(u? — u)!; 
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they are complex for 1 > u > 0 and it follows from (17) that 
(35) A: =1 (u > 0). 


For 1 > u > 0 we therefore have convergence for all g from (0, 2). 
If, on the contrary, u < 0, we obtain 


A, = 1 — 2u+ 2(4 — u)! (u < 0) 


and this exceeds 1. To obtain the value gq; of g between gp and 2 for which 
A, becomes 1, observe that 


N(- 1) = (2-—g? + ug, N(1) = gi — x). 
For u < 0, N,(1) does not vanish while V,(— 1) vanishes for 
(36) gi = 2/(1 + |ul#) (u <0). 


Thus q: lies always between qo and 2, and since for g = q; the product of two 
roots of (13) is (qi — 1)? < 1, we have indeed 


Ag, = 1 (u <0). 


We have therefore for u < 0 convergence if and only if g lies in the interval 
(0, q:). In particular we have here for suitable values of g convergence for all 
negative u, but these values become small for large value of — u. Theorem 3 
is thus completely proved. 

It remains finally to answer the question for what values of g does the in- 
complete relaxation give any improvement at all, that is to say, A, < |u| 
(of course we assume here |u| < 1). 

Now we have, by (14), Ai = |u| and the same is true for g = 1 + u: 


(37) Ai = |14| = Aizu (|u| <1). 


Indeed, we verify immediately that 1 + u = qo according as u is positive 
or negative ; therefore the roots of (13) are complex for g = 1 + u, and (37) 
follows from (17). 

But now it follows from theorem 2 that we have improvement in the case 
of the incomplete relaxation (in contrast to the case g = 1) if and only if g 
lies in the interior of the interval between 1 and 1 + u. Theorem 4 is thus 
proved. 

We discuss finally a special example of a symmetric 2 X 2 matrix in 
which the improvement of the convergence due to the introduction of the 
incomplete relaxation is easily demonstrated explicitly. 


1 
Le a-(2, %). 
We have in this case by (10), (11), and (14) 
(38) Ar = u = 9/25 > 1/3, go = 10/9, Ag, = 1/9. 
The formulas for the cyclic single step iteration with g = 1 are in this case 
(39) xyPHD = — Zx—0)/5, xg) = — 3x, 06+ /5 


and it is readily verified that, in putting —, = (x, x2), we have 


E, = (9/25)” Eo. 
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In taking & = (1, 1) we have then 
(40) \é,| = v2(9/25)’. 


Consider on the other hand the over relaxation with the value of g = 
10/9. Here the components of the approximating vectors are to be computed 
from the equations 


40D) = — 44/9 — 2x2 /3, xt = — x /9 — 2x, 04/3, 


If we put & = (x:, x2) and assume again & = (1,1) = &, we obtain as is 
readily verified 


& = 3-*"1(3 — 24y, 3 + 8y) (vy = 0, 1,---). 


Here we have 
|é,| ~ 8(10)*»9-*/3 (v— o). 


We give in what follows a table of the initial values of |£,| and |é,}. 


v |é lé| 

1 0.5091 0.8780 

2 0.1833 0.2010 

3 0.06598 0.03388 

4 0.02375 0.005048 
5 0.008551 0.0007037 


Although the difference between go and 1 is very small, in fact 1/9, the 
improvement is already observed at £; and becomes more and more pro- 
nounced from there on. 

A. OsTROWSKI 
American University, Washington, D. C. 
University of Basle, Switzerland 


This paper was prepared under a contract of the National Bureau of Standards with the 
American University, Washington, D. C 


1S. P. FRANKEL, “Convergence rates of iterative treatments of partial differential 
equations,” MTAC, v. 4, 1950, p. 65-75. 

2A. M. Ostrowski, On the Linear Iteration Procedures for Symmetric Matrices. NBS 
Report no. 1844, August 1952, p. 23, (68). 


The Accuracy of Numerical Solutions of 
Ordinary Differential Equations 


1. Introduction. The present paper describes a general method by which 
the random and systematic errors may be estimated of numerical solutions of 
any systems of ordinary differential equations. The errors arise from the ac- 
cumulation of rounding-off errors, and from the use of erroneous formulas 
for performing the numerical integrations. The estimation is based on the 
properties of the solutions of the system of equations adjoint to the varia- 
tional equations of the problem, and is applicable to any method of integra- 
tion. 
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The present method was described by the author at a meeting of the 
American Astronomical Society on February 1, 1946 at Columbia Univers- 
ity. Development of the method resulted from a conversation with CHARLES 
B. MorreEy, JR., who explained to the author properties of the solutions of 
adjoint equations with which the author was not then familiar. The pro- 
cedure seems intuitively obvious and straightforward, and it has not been 
published earlier both for this reason and because it was understood that 
Hans RADEMACHER was planning to publish a similar and possibly inde- 
pendent treatment. It now appears, however, that Rademacher! was con- 
cerned with the accuracy of particular methods of integration, while the 
present method is applicable to any integration procedure that may be em- 
ployed. There are still other methods for estimating the accuracy of numeri- 
cal solutions of special types of ordinary differential equations. For example, 
BRouWER? has made special studies of the accuracy of numerical integra- 
tions, by the Crommelin-Cowell method, of the orbital differential equations 
of dynamical astronomy. 

A virtue of the present method is its generality ; but there are alternative 
general methods, possibly just as good, which may not have been published. 
The author understands from conversations with L. H. THomas, that he 
has made use of general procedures not involving the adjoint equations. 
The main justification for publishing a description of the author’s procedure 
is his hope that it may help others to select computational procedures, for 
numerical integrations, that will yield results of desired accuracy. 

2. The Adjoint Equations. Consider the system of first-order differential 
equations 


(1) = Da,a;+5:; ij =1,2,---n 


j=l 
where the ? quantities a;,; and the m quantities }; may vary with the inde- 
pendent variable, ¢. Let \; be a set of variables satisfying the adjoint system 
of equations 
(2) -i= E 4. 
Since 

(d/dt) dea = ota + veh 
(3) = x G;, xi — x @; xij + x bax 

= Dba; 

it follows that 
(4) be (A) \; (A) = x x; (0) A; (0) + f . x bid; dt. 


3. Application. From any system of ordinary differential equations there 
can be derived a set of variational equations of the form (1) where the x;(t)'s 
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are differences, between the exact solution of the original equations corre- 
sponding to the desired initial conditions, and any neighboring exact solution 
not subject to the desired initial conditions. If a single error were made in 
the course of solving the original system of equations by a scheme of stepwise 
integration that was otherwise perfect, the solution would be exact before 
the error was made; for later values of the independent variable the solution 
would still be an exact solution of the original differential equations but 
would correspond to altered initial conditions. The x,’s would all be zero 
before the error, and would grow after it in accordance with the equations of 
the form (1), with 5,’s that were zero for all steps except the one in which 
the error was made. 

If, because of defective methods of calculation or for any other reason 
errors ¢,;(¢) are introduced into the 7’th variable x; at a particular step (¢ — w 
to t, say) in the solution of the original system of differential equations, 
one may consider the ¢’s to have been introduced by },’s in (1) such that 


im f * b(t) dt 


and that are zero outside of the interval (t — w, #). One may consider ap- 
proximately that ; 
bi(t) = (1/w)e:(2) 


throughout the interval and therefore approximately that 
t 
f b(t) (2) dt = e:(t) r,(é). 
t-—w . 


Thus by (4) the resulting final error (at ¢ = A, say) in a particular variable 
(e.g., %1 say) is 


x1(A) = © s(t) Md) 


i=1 


and if errors are introduced at all steps 


(5) x,(A) = > pa e:(t) r,(t) 
all steps i—1 
provided that the \’s are any solution of (2) satisfying the boundary condi- 
tions 
di (A) = 1; A; (A) =0, 7 1. 


4. Truncation and Rounding Errors. Equations (4) and (5) provide 
means for predicting the errors of numerical solutions of systems of ordinary 
differential equations. Rough solutions of the equations (2), based on a 
rough solution of the original equations, are in practice adequate. Rounding- 
off errors of the usual hand-made variety are drawn from populations whose 
means are zero, and whose individuals are uniformly distributed from —1/2 
to 1/2 in units of the last digit. Their variance is thus 1/12 in such units. 
The resulting variance of a final value like x, (A) is given by 


(6) 075, (A) =o; >> Ar (2) + o? > d2?(t) +---+6,2 >, dw2(t) 
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where a? is the variance of the rounding-off errors introduced into the vari- 
able x; at any step, and where the sums are taken over all steps. The preced- 
ing equation is general; if the rounding-off errors are not of the usual hand- 
made sort it is still valid. If the rounding-off errors do not come from popula- 
tions with zero means, then a bias, or systematic error is introduced whose 
final value has the population mean 


H(A) = DT MO) Alt 


where M(t) is the population mean of the error ¢;(#), conceivably a function 
of ¢. 

Besides rounding-off errors, ‘‘truncation’’ errors are introduced by the 
circumstance that the formulas employed for integrations are erroneous. 
Whatever the formulas are, and however they are employed, iterated or not, 
any particular method of integration applied to a particular system of differ- 
ential equations always corresponds to the exact solution of a system of 
difference equations rather than of the original differential equations. The 
particular method of integration thus corresponds to an exact solution of a 
system of differential equations somewhat different from the original differ- 
ential equations. It is always possible to evaluate, approximately, the differ- 
ences between the original differential equations and those that the scheme 
is exactly solving, then to find the appropriate 0,’s in the variational equa- 
tions of the form (1), and finally to apply (4) to predict the final errors, thus 


(7) ni(A) = fo Eb asl at 


in which, as usual, the \’s must be chosen to satisfy the boundary conditions 


Alternatively, one can find appropriate truncation errors ¢;(#) and then 
apply equation (5). 

For planning purposes, no great accuracy in the calculations of accuracy 
is necessary, and no great accuracy should be sought. Even rough calcula- 
tions are expected to suffice to decide how many digits, what size of steps, 
and what scheme of integration to employ. 

5. Example. An example, suggested by WERNER LEUTERT, will be given 
of the use of the preceding method by an application to the non-linear 
differential equation of the first order 


(8) y = (3/2) ty"* 


which can be integrated analytically but which will be treated as though it 
can not be. Suppose that one wishes to integrate this equation numerically, 
starting with the value 


y(i) = 1, 


as far as y(5); and that one wishes the value y(5) to be accurate ‘‘to the 
third place” of decimals. One wishes to use the scheme of integration defined 
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by the approximate formula 
Vv v? " 
(9) vy = [1-5-7 foo 
in which the differences are backward, or ascending, and in which w is the 
length of a step. One wishes to determine w, and the number of decimals to 


retain in the calculations. 
The variational equation corresponding to (1) is 


(10) b= 5ty Me 

and the adjoint equation corresponding to (2) is thus 

(11) } = Sty. 

A rough integration of (8) must first be accomplished, and then a rough inte- 
gration of (11) to find A. Such integrations have been accomplished by the 


aid of a ten-inch slide rule, and the results appear in the first four columns 
of the following table: 


t y 2.27A r ADty 
1 1 1 44 25 
2 2.80 1.44 .63 .06 
3 5.15 1.76 .78 .03 
4 7.95 2.03 89 01 
5 11.15 2.27 1.00 01 


No attempt has been made to obtain results correct to the second place of 
decimals, although two places were retained. The integration for \ with the 
starting value unity led to a value 2.27 at ¢ = 5; the fourth column contains 
d adjusted to have the value unity at ¢ = 5. 

Consideration of the difference formula (9) shows that it corresponds 
substantially to the differential equation 


(8’) Dy — (1/24) w* Dty + (11/720) wt Dy +--+ = (3/2) ty™* 
instead of to equation (8), so that approximately 


b(t) = (1/24) w® Dt y 
™ e(t) = (1/24) wt Dt y; 


these results could have been obtained directly from the term of lowest order 
that has been omitted from the right-hand member of equation (9). By 
equation (5) or equation (7) the truncation error at ¢ = 5 is 


x(5) = (w*/24) f D* y(t) (2) dle 
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Values of D*y are obtained from equation (8); values of D*y(#) A(t) are 
tabulated above; a quadrature furnishes the result 

(12) x(5) = w*/120. 


It is noticed that a unit error in wy at any stage introduces a unit error in 
y when the scheme of integration is (9). Hence by equation (6) the variance 
of y(5) arising from rounding errors is 


(1/12) E2X@) 


t=1 


o7y(5) 


Il 


(1/12 w) f d2(t) dt 
1/5 w 


by quadratures, very nearly, in units corresponding to the last place retained 
in wy. 

To obtain a value of y(5) accurate ‘‘to the third place” of decimals, one 
equates the right-hand member of equation (12) to 0.0005 and solves for w. 
It is found that w is .39. One can adopt a value w = .4, if one will tolerate a 
bias error of 0.00053 in y(5); this is tolerated, and the value w = .4 is 
adopted. A number fairly rounded to the nearest 0.001 has a rounding error 
whose variance is 1/12 in units of the sixth place. It is therefore reasonable 
to require that the variance of y(5) from the accumulation of rounding 
errors should be smaller than 1/12 in the sixth place. If only three decimals 
were retained in the values of wy then the variance of y(5) would be 1/5w or 
1/2 in the sixth place, which is too large to be acceptable. With four decimals, 
the variance is 1/2 in the eighth place, or 1/200 in the sixth, which is better 
than is needed. Therefore four places of decimals should be retained in 
values of wy. 

The definitive integration of equation (8) by the scheme (9) was next 
accomplished, with steps of 0.4 and with four places of decimals. The value 
obtained was 

y(5) = 11.1810. 


The correct value of y(5), obtained by analytical solution of (8), is 


(5) = 5 
= 11.18034------ 


showing that the error of y(5) obtained by the numerical integration is 
0.00066, and that (5) is substantially as accurate as was desired and as was 
predicted. 

THEODORE E. STERNE 
Ballistic Research Laboratories 
Aberdeen Proving Ground 
Maryland 


1 See, for instance, HANS RADEMACHER, “On the accumulation of errors in numerical 
integration on the Eniac,” lecture 19 of Theory and Techniques for Design of Electronic 
Digital Com s, v. 2, Moore School of Electrical Engineering, University of Penn., 1947. 

2 Dirk Brouwer, ‘‘On the accumulation of errors in numerical integration,” A stronomi- 
cal Jn., v. 46, p. 149, 1937. 





THE ZEROS OF THE PARTIAL SUMS OF é? 165 
are The Zeros of the Partial Sums of e’ 
The location of the zeros of certain entire functions has received con- 
siderable attention in the literature, and it was suggested by R. S. VARGA 
‘in that a table of the zeros of certain truncated power series might be of interest. 
nce F . nm gh 
Accordingly, the zeros of the truncated exponential series S,(z) = > ki 
> FR! 
have been computed for values of m up to 23. They appear in the accompany- 
ing table with the complex zeros ordered as to modulus and with the real 
zero printed last. 
Table of Zeros of S,(z) 
n=2 n= 11 
— 1.0000 0000 0000 + 1.0000 0000 0000: — 3.6641 5160 2429 + 1.5570 4385 7798: 
1ed — 2.9170 5082 0338 + 3.0563 7182 98311 
n= 3 — 1.5814 pon a + — maak 7709 
a 0.5460 781 + 5.524! 6051i 
one ern. ee 4.0692 9095 3063 + 6.0474 9237 98891 
w. ’ —3.9054 5175 7616 
ea n=4 
- is — 1.7294 4423 1067 + oo ese saan 2 = 12 
—0.2 2 « i , 
ror aatsdirns ind en * 4.1356 0823 9264 + 0.7755 4204 9110: 
ble aes ~saees 9710 2446 + ase ores 9247i 
i . — 2.7579 8923 1191 4 3. 23534 
ing — 1.6495 0283 1735 + 1.6939 3340 4349i = __ 12491 2514 3358 4 5.0495 5510 7284: 
als 0.2398 0639 3753 + 3.1283 3502 5970: 1.0534 2363 9656 + 6.0594 9143 3864: 
) or — 2.1806 0712 4035 4.7781 9607 6918 + 6.4511 7633 7446: 
als, ~~ y 
ter — 2.3618 1018 0482 + 0.8383 $027 79174 n= 13 
in — 1.4418 0139 0549 + 2.4345 2268 1808: —4.2712 4352 2040 + 1.5348 5553 2416i 
0.8036 1157 1031 + 3.6977 0175 3629% — 3.6448 0690 0313 + 3.0273 4405 4885; 
ext ; —2.5489 2149 0908 + 4.4259 0791 2207: 
1 _— —0.8813 4146 1503 + 5.6544 1733 8744: 
tue — 2.3798 8388 3168 + 1.6289 9897 6372: 1.5843 14949756 + 6.5740 0727 9114: 
— 1.1472 0068 9937 + 3.1240 3923 8058: 5.4997 0440 2346 + 6.8391 5923 4366: 
1.4065 8592 8087 + 4.2250 6684 49497 — 4.4754 1195 4676 
— 2.7590 0270 9962 
n= 8 n= 14 
— 2.9645 9950 5160 + 0.8088 7832 7313: —4.7125 8682 7652 + 0.7651 0266 4661: 
— 2.2864 2928 4171 + 2.3777 1166 77931 —4.3307 0823 3773 + 2.2772 3593 3932: 
—0.7887 9362 0387 + 3.7718 1078 3950: — 3.5438 3446 4412 + 3.7319 2259 2644; 
. 2.0398 2240 9719 + 4.7186 1488 3923: — 2.2976 9841 1869 + 5.0783 9242 5805: 
i is —0.4831 5856 3940 + 6.2391 4928 3405: 
vas n=9 2.1356 7746 7731 + 7.0705 6793 5765i 
— 3.0386 4807 2936 + 1.5868 0119 5758i 6.2525 COGS SPIT te 7.2151 4805 BODES 
— 2.1108 3980 0302 + 3.0899 1092 8725: 
— 0.3810 6984 5663 + 4.3846 4453 31454 15 
2.6973 3346 1536 + 5.1841 6206 26497 oe 
— 3.3335 5148 5269 —4,8669 6227.5703 + 1.5176 2913 7005: 
—4.3271 6518 2166 + 3.0027 7099 6726: 
Fy Spee Lan oh te 
ical — 3.5538 7599 3928 + 0.7894 2208 2895: — 2.0103 3 + 5. 49 17373 
_~s — 3.0155 3577 0425 + 2.3352 2385 7750: —0.0585 5212 5116 + 6.8056 2430 67914 
947. — 1.8716 6001 0419 + 3.7701 9023 1409 2.7050 4956 8386 + 7.5509 4048 4655: 
ym1- 0.0662 0154 6301 + 4.9676 7937 0404: 6.9747 8093 2988 + 7.5745 6159 4666: 
3.3748 7022 8472 + 5.6260 2017 96981 — 5.0438 8707 2612 
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Table of Zeros of S,(z)—Continued 


— 5.2863 1780 3783 + 0.7569 4362 0482: — 6.4273 0264 0861 + 0.7449 7139 42103 
—4,9527 8138 2083 + 2.2566 6238 4481i — 6.1611 0097 0908 + 2.2255 2308 T4174 
—4.2704 2428 9874 + 3.7119 3355 5596i —5.6209 8880 2007 + 3.6770 7242 0762: 
—3.2047 4966 4526 + 5.0858 8483 8575i —4.7903 2735 9920 + 5.0773 1588 8906: 
—1.6915 4717 7729 + 6.3274 3910 8059i — 3.6406 2228 5764 + 6.3986 7391 19691 
0.3892 9335 7090 + 7.3554 4605 9723: — 2.1255 4693 1251 + 7.6040 5463 9680: 
3.2904 2476 4254 + 8.0166 1911 0269: — 0.1684 0640 7437 + 8.6388 4961 07321 
7.7261 0219 6653 + 7.9245 9187 5073: aaa _ = + 9.4133 5826 0670: 
: 7 + 9.7554 8801 5501: 
10.8045 8424 5908 + 9.2291 9790 46771 

n= 17 n= 21 
—6.6167 8934 1171 + 1.4830 8180 56311 
—5.4551 0328 9602 + 1.5038 3976 7334: —6.2346 1169 3725 + 2.9488 4393 55021 
— 4.9806 7273 2968 + 2.9819 4653 1042: —5.5863 1103 9618 + 4.3785 4166 8723: 
— 4.1680 2106 9422 + 4.4053 7601 56391 = 4.6531 2333 7897 + 5.7502 8769 4696: 
— 2.9788 2514 8596 + 5.7375 9875 5455i —3,4046 1507 1622 + 7.0364 9649 65327 
— 1.3451 2663 7403 + 6.9268 7649 60751 — 1.7924 6476 9782 + 8.1996 0590 2952: 
0.8577 8157 7290 + 7.8899 9871 2830: 0.2624 14640116 + 9.1839 0151 88837 
3.8901 4238 0621 + 8.4688 8076 9632 2.8999 6466 8596 + 9.8974 9663 3020: 
8.4854 1859 0153 + 8.2642 5429 21761 6.4075 0945 9944 + 10.1637 8334 0590: 
—5.6111 8734 0141 11.5895 7008 7260 + 9.5350 4977 9337: 

— 6.7430 9089 3669 

18 n = 22 
us —6.9955 1940 2365 + 0.7404 3625 9274i 
—5.8576 9828 4668 + 0.7503 7782 8924: —6.7537 1576 2515 + 2.2134 4348 11017 
—5.5616 1575 6798 + 2.2397 2181 6511: — 6.2643 4318 6033 + 3.6622 7895 1044: 
—4.9588 10200174 + 3.6935 9122 5275: —5.5151 2445 5329 + 5.0687 9151 8941: 
—4.0258 8765 2371 + 5.0838 2295 5316: —4.4855 3901 9504 + 6.4113 1616 40467 
—2.7214 0644 5610 + 6.3738 9942 0060: — 3.1434 8080 8797 + 7.6621 3795 8098: 
—0.9741 5991 2999 + 7.5112 3498 4979; — 1.4388 3815 7314 + 8.7831 4914 18247 
1.3447 6355 9230 + 8.4104 8615 4221: 0.7096 9887 5451 + 9.7174 9498 5677: 
4.5028 0973 4285 + 8.9088 2705 4271: 3.4453 7712 5688 + 10.3711 1273 05127 
12.3797 8501 1070 + 9.8339 1911 3179% 

n = 23 
n= 19 —7.1926 8590 7451 + 1.4750 5919 5082i 
—6.0379 0688 9211 + 1.4925 34019264;  — O-B483 LALI BOGS & 2.9355 2075 2708 
—4,8936 3517 6394 + 4.3919 0787 8943: — 42905 2911 0822 + 7.0608 9517 9930: 
— 3.8487 8978 9420 + 5.7480 1406 1058: ' ; . 
34360 0924 7000 . 6.9957 S579 15477 — 2.8595 2302 8914 + 8.2762 1397 0162i 
a chineits 5 pense. t — 1.0664 4595 7935 + 9.3553 6015 89097 
—0.581 + 8.0815 7686 4150: 1.1720 9815 7227 + 10.2403 1777 0398: 
1.8484 3722 1795 + 8.9179 6283 4128: 4.0025 2534 8326 + 10.8348 6485 1671: 
5.1272 4541 2251 + 9.3374 1621 1924: 7.7243 3865 4741 + 10.9536 0414 4523: 
10.0252 3949 6630 + 8.9158 4881 41227 13.1748 6483 8907 + 10.1262 6417 8727: 

—6.1775 3407 8278 + —7.3079 8221 4646 


The work was done on the Harvard Mark IV Calculator! using the 
quadratic factor method? together with a routine which supplied initial ap- 
proximations chosen on a rectangular mesh in the complex plane. Mark IV 
operates with a fixed decimal point and carries sixteen decimal digits, but 
to reduce round-off errors a floating point routine was used. As a final check, 
the sum and product of the roots of each polynomial were compared with the 
appropriate coefficients of S,(z), and agreement to twelve significant digits 
was obtained for all m < 21. 
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The plot of the zeros in Fig. 1 exhibits the regularity of the family of 
curves joining the m zeros of a given S,(z). The broken curves which join 
zeros of the same rank (when ordered as to modulus for each nm) have immedi- 
ate application in that the zeros of the partial sum S,(z) may be located ap- 
proximately from a knowledge of the zeros of the partial sums of lower order. 
This could be used to provide good first approximations in extending the 
table. 





Fic. 1. 


It is interesting to consider the present numerical results in the light of 
the following known properties of S,(z) : 


(1) Lim = = ; _ ~ where r, is the number of zeros of S,(z) lying in 
the right half-plane.* 

(2) Sen(z) and Seas1(z) each have 2m complex zeros.‘ 

(3) The semi-infinite strip® |y| < ¥6, x > 0, contains no zeros of any 
S,(z). Figure 1 suggests that much larger zero-free regions exist. 

(4) Every zero of S,(z) satisfies the inequality® 


>|s| >> 
n 2 a: 


(5) The zero of smallest modulus of S2,,:(z) is real and negative. An 
unpublished proof has been given by D. J. NEWMAN. 

(6) The convex polygon formed by the zeros of S,(z) encloses all the 
zeros of the partial sums of lower order. This is based on the fact 
that S,’(z) = S,-:(z) and on the following theorem :? Any convex 
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polygon which contains all the zeros of a polynomial p(z) also con- 
tains all the zeros of the derivative p’(z). 
K. E. Iverson 
Harvard University 
Cambridge, Mass. 


14 ae of the Mark IV Calculator. Computation Laboratory of Harvard Uni- 
= -- v. 28 (In preparation). 

W.E. Mune, Numerical Calculus. Princeton 1949, p. 53; or D. R. HARTREE, Numerical 
Ausipal Oxford 1952, p. 205. 


3G. SzEG6, “(ber eine Eigenschaft der Exponentialreihe,’’ Berliner Mathematischen 
Gesellschaft, Sitsungsberichte, 1924, p. 50-64. 


= 3. BERGHUIS, “Truncated power series,” Mathematical Centre, Amsterdam, Report 
R173, 1952. 


5 R. S. VarGa,. “Semi-infinite and infinite strips free of zeros,’’ Universita e Politecnico 
di Torino, Seminario Matematico, Rendiconti, v. 11, 1951-1952, p. 289. 


6K. S. K. Iyvencar, “A note on the zeros of » = Mathematics Student, v. 6, 1938, 


p. 77. 
7M. MarRpDEN, The Geometry of the —_ a a Polynomial in a Complex Variable. Ameri- 
can Mathematical Society. New York, 1949 
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1094[B ].—Marcuant CALCULATING MACHINE Co. Table No. 80 of Factors 
for 8-Place Square Roots, 1951, 4 p.; Table No. 81 --- for 6-Place Roots, 
1952,4 p.;and Table No. 82 --- for 5-Place Roots, 1952, 2 p. Publications 
of Marchant Calculators, Inc., Oakland, California. 21.5 * 28 cm. 


The tables for 5- and 6-place roots resemble a former one for 5-place roots 
of the same publisher (MTAC, v. 1, p. 356). As in the former table the root 
results from adding a tabular number to the number WN of which root is de- 
sired, and dividing this sum by an adjacent tabular number. The table for 
8-place roots requires two divisions but without need of intermediate copy- 
ing. Mathematically the table for 8-place roots is equivalent to a similar 
two-division table reported in MTAC, v. 5, p. 180. 

Doubtless the most welcome of these tables will be the one for 6-place 
roots because this number of places is specified by the majority of work 
sheets for general surveying and military uses. There has not been available 
a means of obtaining 6-place roots by the use of tables of this type and a 
single division. 

Significant improvement has been made in tables Nos. 81 and 82 as com- 
pared with the previous 5-place table (1) by eliminating need of determining 
which of two arguments is nearer to N, (2) by reducing the number of tabu- 
lated divisors—resulting from using a single column of divisors for the entire 
range of NV from 1 to 100 and from altering the range of error from (0 to +5) 
to (—5 to +5) in units of last place, (3) by using the leeway afforded by 
integer arguments to reduce the number of digits of the divisors to as few as 
possible. 

The method of selecting argument intervals and the error expression of 
tables of the single-division type are described in Willers, Practical Analysis, 
Dover edition, Art. 6-21, 1948 (MTAC, v. 3, p. 493). Obvious alteration of 
the error expression is necessary when applying them to these tables because 
of change of the range of error. 
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The tables are suitable for use with any of the usual types of desk cal- 
culating machines, though the manipulative description accompanying table 
80 refers to the Marchant machine. The tables were prepared under the 
direction of H. T. Avery who asks that aid rendered by H. J. REyNoLpDs be 
acknowledged. 

T. W. Smmpson 
66 Alvarado Road 
Berkeley, California 


1095[B].—Victor THeBAuLT, Les Récréations Mathématiques (Parmi les 
nombres curieux). Paris, Gauthier-Villars, 1952, vi, 299 p., 15.5 XK 23.5 
cm. 


P. 12: Table of the 87 squares formed by the digits 0, 1, 2, 3, 4, 5, 6, 
7,8,9, each once and only once. 
P. 89-119: Tables of squares of integers 1(1)1000 for bases 2(1)9, 11, 12. 
P. 226-230: The tables of squares with (a) 3, (b) 5, (c) 6, (d) 7, (e) 8, 
(f) 9 digits occurring once and only once. 
R. C. ARCHIBALD 
Brown University 
Providence, R. I. 


1096[C,K ].—M. S. BartLett, “The statistical significance of odd bits of 
information,” Biometrika, v. 39, 1952, p. 228-237. 


Table I (p. 230) lists to 4D for » = 0(.01)1 the functions — In p; 
— pln p; pln? p; — pln p — (1 — p) In(1 — p); and p(1 — p) In* (p/(1 — p)). 
An example is given using this table to compute an “information” function. 


W. J. Dixon 
University of Oregon 


1097[C ].—ALEXANDER JOHN THOMPSON, Logarithmetica Britannica, being 
a Standard Table of Logarithms to Twenty Decimal Places. Part II, Num- 
bers 20,000 to 30,000 together with General Introduction. Cambridge Uni- 
versity Press, 1952. vi, [102], xcviii, [xi] p., 21.6 & 27.7 cm., 45 shillings; 
New York price $8.50. 


This is the ninth and final part of the very notable work which Dr. 
Thompson started 30 years ago, as a publication for KARL PEARSON’s Tracts 
for Computers, Department of Statistics of the University of London, to 
commemorate the tercentenary of HENRY BricGs’ publication of Arithmetica 
Logarithmica, 1624. This work of Briggs (1561-1631) contained log N to 14D 
for N = 1(1)20000, 90000(1)101000, and the square roots of the integers 
[1(1)200; 11D], with first differences in each case; (see MTAC, v. 1, p. 170; 
v. 2, p., 94,196). On pages Ixxviii-Ixxxiii of this part II is a list of errors of 
more than a unit in the values of log N in this Briggsian table. Briggs had 
practically completed the computation of log N for N = 20001(1)90000 
when ApRIAN VLAcQ published Arithmetica Logarithmica (1628) for log N, 
N = [1(1)100000; 10D], and characterized his volume as the second edition 
of the work by Briggs. The completion by Vlacq of a 10D table from the 
earlier 14D tables of Briggs was a comparatively simple matter. But his 
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action in rushing into print with such a publication can scarcely be termed 
other than reprehensible. 

The first part of Dr. Thompson’s work, N = 90,000—-100,000, appeared 
in 1924, and this was followed by other parts in 1927, 1928, 1931, 1933, 1934, 
1935, 1937. Thus 15 years elapsed after the appearance of 8 parts, before the 
final part, completing the table for N = 10000(1)100000, was published. 
For brief reviews of parts VI and VII see MTAC, RMT 23, 42, 65, v. 1, 
p. 4, 5, 7, that is, Scripta Mathematica, v. 2, 1934, p., 196; v. 3, 1935, p., 192; 
v. 4, 1936, p., 201. 

Since 100 logarithms are displayed on each page, there are in the whole 
work 900 pages occupied by the fundamental table. To introduce this table 
there are in the present part 108 pages of preliminary matter. In other parts 
there was further preliminary matter to which we shall presently refer. 
Hence in this part are the necessary title-pages and indices for binding the 
whole work into two volumes. Of the useful factoring Table F (logarithms of 
1+ N/107, 1 + N/10", and 1 + N/10"%, N = [0(1)1000; 21D ]}) there are 
two copies, one for each volume. Reprints of the complete work in two bound 
volumes are now available at the publisher. For completing sets, copies of 
the first eight published parts are available at 21 shillings each. 

In part I (1934) is a two-page facsimile of the original will of HENRY 
BriGcs, the first page in his own handwriting; part III (1937) has four pages 
with facsimiles of five letters of Briggs, and a Note by Dr. THompPson; in 
part IV (1928) there is a facsimile of the title-page of the work in which 
Briggs’s Treatise on the Northwest Passage to the South Sea was published 
in 1622; two pages of errors in Arithmetica Logarithmica, 1624, are also given ; 
part V (1931) has a facsimile of the title-page of Arithmetica Logarithmica ; 
and part VI (1933) has a facsimile of a Briggs’ letter to JoHN PELL. 

The Mirifict Logarithmorum Canonis Descriptio was first published in 
1614 by Joun Napier (1550-1617). The second edition in 1619, prepared 
after Napier’s death by his son ROBERT, and Briggs, contained an edited 
earlier work of Napier : Mirifict Logarithmorum Constructio, Canonts in which 
Napier called logarithms ‘‘artificial numbers.’’ See my article, ‘“‘Napier’s 
Descriptio and Constructio,’’ Amer. Math. Soc., Bull., v. 22, 1916, p. 182-187. 
Part VII (1935) has four pages (three of 6 facsimile pages) illustrating the 
relation of Henry Briggs to Napier’s Constructio Canonis with a Note by 
Kart Pearson. It is by no means clear how Pearson reasons that the 
word “logarithm” used by Napier in his work of 1614 “must have been 
added by Briggs or Robert Napier,”” because Napier had used a different 
term in his Constructio. 

Parts VIII (1927) and IX (1924) contained a number of pages of intro- 
ductory matter by Thompson and in the latter a ‘‘Prefatory Note” by Karl 
Pearson, and also a facsimile of the title and specimen pages of the excessively 
rare 1617 tract (16 p.): Logarithmorum Chilias Prima of HENRY BrIGGs. 
This is the earliest publication of logarithms to the base 10; it exhibits the 
logarithms to 14D of numbers 1(1)1000. Table F was given in part VIII. 

Dr. THompson’s splendid Introduction to Logarithmetica Britannica in 
part II includes a discussion of: (a) Interpolation with examples (p. xv— 
xxviii); (b) Method of Construction (p. xxviii-lxii)—here are 12 special 
tables and a plate illustrating the unique four-bank integrating and differ- 
encing machine which Thompson built for carrying through the calculations 
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of this work. He personally set up the type for all table entries by means of 
a monotype keyboard. This system of type setting involves the use of two 
entirely separate machines, a keyboard and a typecaster. Thompson used 
the keyboard to punch holes in a continuous ribbon of specially prepared 
paper, which then went to the typecaster. 

There is a section (p. lvi—lxii) on construction of a table of Anti-logarithms. 
Shortly after the publication of the table of logarithms had begun, Thompson 
was urged to undertake a companion volume of anti-logarithms on the same 
scale. Table 8 (p. Iviii) is a single page of such a table log N = 0.00(.01)1.00, 
corresponding N’s being given to 28D. 

“Printing and proof-reading’’ and “Acknowledgments” on p. Ixii-Ixv 
are followed by ‘“‘References’’ (p. lxv—lxvi) listing works mentioned in the 
Introduction, and a few other titles. 

In the manuscript for the main table the calculations were made to about 
24D, the last digit having a possible error of three or four units. In starting 
computations it was desirable to have a considerable number of values of 
log N to a large number of decimal places. For one thing Dr. THompson 
checked by 63D computation the accuracy of the 61D table of log N, N = 
1(1)100 and primes to 1097; also 999990(1)1000110, in ABRAHAM SHARP’s 
Geometry Improv'd (1717). The figures in this table are arranged in five-figure 
groups except the twelfth and last group, which consists of six figures. The 
following errors in mantissas were found (the error for N = 751 was new): 


N Group For Read 
103 9 33496 23496 
227 12 494656 495656 
751 12 287788 287771 
839 12 539741 538741 
1009 12 382385 382285 


No other errors were found except for unit errors in the 61st decimal place. 
For N = 127, 149, 293, the final digits should each be increased by unity. 

Thompson had only the CALLET (1795) reprint by Sharp, and found that 
Callet had an error for N = 1097; but here Sharp was correct. In the PETERS 
& STEIN reprint (1922) all of Callet’s errors persist. 

Other extensive log N tables are: (a) A. GRIMPEN’s 84D table (1922) of 
the prime numbers up to 113; and (b) H. M. Parkaurst’s 100D table (1876) 
for 96 values of N < 109. See MTAC, v. 1, p. 20, 58-59, 121-122. 

Appendix (i) of part II contains the first English translation of the 
earliest published biography of Briccs, in THoMAs Smitu’s Vitae quorundam 
eruditissimorum et illustrium virorum, 1707. There is no known portrait of 
Briggs (see MTAC, v. 2, p. 287; v. 3, p. 67). 

Appendix (ii) lists the errors in Arithmetica Logarithmica (1624), referred 
to above. 

Then follow: Table F (Ixxxv-xciv); Table G, antilogarithms of loga- 
rithms [0000000 (1)0000450 ; 21D], p. xcv—xcvii; Table H, Short Tables and 
Constants, p. xcviii; and finally (p. xcix-cv) log N, N = [1(1)1000; 21D]. 

The nine parts of this work have been published in Tracts for Computers 
prepared in the Department of Statistics, University of London, in the fol- 
lowing numbers: XIX, XXII, XXI, XVI, XVII, XVIII, XX, XIV, XI. 
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Dr. THompson tells us that during the past 30 years no error has been 
reported for any of nearly 3000000 figures in the eight parts published up to 
1937. We tender the author our heartiest felicitations on his truly monu- 
mental personal contributions achieved in producing this great work. 

R. C. ARCHIBALD 
Brown University 
Providence, R. I. 


1098[F ].—H. J. A. Duparc, C. G. LEKKERKERKER & W. PEREMANS, 
Reduced sequences of integers and pseudo-random numbers, Math. Cen- 
trum, Amsterdam, Report ZW1952-002, 15 mimeographed leaves. 


This study of the length of the period of a geometric progression or a 
Fibonacci sequence when the terms are reduced modulo m, contains the 
following small tables. 

P. 3 contains a table of the factors of 10" — 1, Euler’s totient function of 
these factors and their least common multiple for » = 1(1)10. The corre- 
sponding information for 10" + 1 is given for m = 1(1)7. 

P. 4 gives the same information for 2" + 1 for m = 1(1)16, 29, 30. 

P. 15 gives the rank of apparition of p in the Fibonacci sequence 0, 1, 1, 
2, 3,5, --+, for each prime less than 433. This table was taken from the 
table of JARDEN.! 

D. H. L. 


1D. JARDEN, ‘Table of the ranks of apparition in Fibonacci’s sequence,” Riveon Lemate- 
matika, v. 1, no. 3, 1946, p. 54 [MTAC, v. 2, p. 343]. 


1099[F. }—Karit GoLpBerG, “‘A table of Wilson quotients and the third 
Wilson prime,” London Math. Soc., Jn., v. 28, 1953, p. 252-256. 


The Wilson quotients W, are defined as the non-negative residues modulo 
p of [(p + 1)! + 1]/p, where p is prime and the Wilson primes are solutions 
of the equation W, = 0. The present table gives the values of W, for all 
primes less than 10000, and shows 563 to be the third Wilson prime. 

N. G. W. H. BEEGER’s table! of Wilson Quotients extended to all primes 
less than 300 and the previously known Wilson primes were 5 and 13. 

GOLDBERG states that six months after the completion of this table 
DonaLD WaLL [UMT 150, MTAC, v. 6, p. 238] computed a table of Wilson 
quotients for all primes less than 5000, and that his table checked with this 
one in every case. 

R. C. ARCHIBALD 

Brown University 
Providence, R. I. 


1N. G. W. H. BrEEceErR, “On the congruence (p — 1)! =— 1 (mod #°),”” Messenger 
Math., v. 49, 1920, p. 177-178. 


1100[F ]..—F. GRUENBERGER, Table of Prime Numbers from 2 to 406253. 
Numerical Analysis Laboratory, Univ. of Wisconsin, Madison, 1953, 1 
“microcard,” 7.6 X 12.4 cm. Price, 25 cents. 


This little card contains 34320 primes, each prime fully spelled out, as 
photographed from 65 sheets, calculated and printed with an IBM Card 
Programmed Calculator. 
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The condensation achieved is remarkable. Twenty cards of this size 
would accommodate the 665000 primes in LEHMER'’s list! covering the first 
10 million numbers. 

The card requires a hand lens of power > 5. The time required to enter 
the table is nearly the same as for a standard table, the time spent in locating 
the appropriate region of the card being comparable with that spent in 
turning pages. 

The list was compared at 520 places (every 66th prime) with that of 
Lehmer. 

Mm. &. 


1D. N. Lesmer, List of Prime Numbers from 1 to 10006721, Washington, 1914. 


1101[F ].—Gruseprre PaLtamMA & L. Potetti, “Tavola dei numeri primi 
dell’intervallo 12 012 000 — 12 072 060,” Unione Matematica Italiana, 
Bollettino, s. 3, v. 8, 1953, p. 52-58. 

With references to the D. N. LEHMER list of primes up to 10 006 721 
(1914) ; and to the list of N. G. W. H. Beecer, L. PoLetti, A. GLopen & 
R. J. Porter, from 10006 741-10 999 997 (1951, MTAC, v. 6, p. 81-82); 
the authors add 3684 primes before 12 072 060. 


R. C. ARCHIBALD 
Brown University 
Providence, R. I. 


1102[F ].—G. Ricci, “‘La differenza di numeri primi consecutivi,” Univ., 
Politec., Torino, Rend. Sem. Mat., v. 11, 1952, p. 149-200. 
This topical history contains the following tables and graphs. 
Let p; = 2, pe = 3, ps = 5, --- be the primes in increasing order and 
let the difference between consecutive primes be denoted by 


d(n) = Pati — Pn. 


Let A»(x) denote the number of primes » < x for which p + 2h is also a 
prime. Let B(x) denote the number of those primes p, < x for which 
d(n) < 2h. 

The main table (p. 152-155) gives the values of 


for m = 1(1)170. The function d(m) is graphed on p. 157 and compared with 
In n. The functions 
Aos(x), (h = 1(1)7) 
and the function 
A(x) =d(m) (bn £ x < Poy) 
are graphed on p. 158-9 for 0 < x < 1020. 
WEsTERN’S' table of those primes p, < 10’ whose difference d(m) exceeds 


that of all smaller primes is reprinted. 
D. H. L. 


1 A. E. WESTERN, “Note on the magnitude of the difference between successive primes,” 
London Math. Soc., Jn., v. 9, 1934, p. 276-278. 
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1103[1,L ].—NBS Tables of Chebyshev Polynomials S,(x) and C,(x). Applied 
Math. Ser. No. 9. U.S. Government Printing Office, Washington, D. C., 
1952. xxix + 161 p., 20.5 X 27 cm. $1.75. 


“The Chebyshev polynomials are of use in many mathematical investi- 
gations. Although direct numerical tabulation is fairly easy to avoid—for 
example, by double or multiple use of ordinary trigonometrical tables—the 
present tables are welcome because they will remove the necessity for these 
roundabout methods, which are often irritating” (from the Foreword by 
J. C. P. MILLER). 

The polynomials tabulated in this volume are 


C,(x) = 2 cosn0@ 


sin (n + 1) 0 
Sa(x) = sin 6 
where 
x = 2cos@. 


Other notations are 
i 
T(x) = 3 C.(2x), T.*(x) = 5 CnC — 2) 
Un(x) = S,(2x), Un* (x) = Sa(4x — 2). 


Explicit expressions are given for C,(x), Sn(x), Tn(x), Un(x) for n = 0 
(1) 12; and for T,,*(x) for = 0 (1) 20. The expansions of x* in the T,(x), 
and in the T,,*(x), polynomials are also given for k = 0 (1) 12. 

The principal tables give 12D values of S,(x) and C,(x) for m = 2 (1) 12, 
x = 0 (.001) 2. 

The Introduction, by CorNELIus LANczos, describes the basic properties 
of Chebyshev polynomials, their applications to expansions, curve fitting, 
solution of linear differential equations with rational coefficients; it gives an 
account of the computation of the tables and instructions for their use, and 
also a list of references. 

The computations were carried out under the technical direction of 
A. N. Lowan. 

A preliminary MS of these tables was described by Lowan in MTAC, 
v. 1, p. 125 (UMT 11). Other tables of Chebyshev polynomials are referred 
to in MTAC, v. 1, p. 385 (RMT 185), v. 2, p. 256 (RMT 371), v. 2, p. 262- 
263 (RMT 381), v. 2, p. 266 (RMT 383), v. 3, p. 97 (RMT 495), v. 3, p. 
119 (MTE 126), v. 3, p. 120-121 (UMT 68). 

“The tabulation of these polynomials—easy to calculate, and easy to 
sidetrack at the cost of some inconvenience—is long overdue. The Computa- 
tion Laboratory staff is to be congratulated on the removal, at last, of this 
source of inconvenience, and in doing so, on the addition of yet one more 
table to its already magnificent series’ (from the Foreword). 

A. E. 


1104[1].—H. E. Sauzer, “Formulas for numerical differentiation in the 
complex plane,” Jn. Math. Phys., v. 31, 1952, p. 155-169. 


In an earlier paper! the author gave coefficients for numerical integration 
in the complex plane of polynomials of degree no higher than eight, based on 
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configurations of the grid points chosen so as to be convenient for initiating 
a computation and “‘as close together as possibie.’’ Thus let the Lagrange 
polynomial of degree m — 1 be defined by 


(1) f(s) = ¥ Pele) fle) 
k=1 


where f(z,) is the given value of f(z) at z, = x, + ty,. The points are chosen 
over a square grid, namely 
2; = 20 + jh 


where j is a complex integer. As an example, for a polynomial of degree four 
the selected grid-points are 


20, 20 a h, Zo + 2h, Zo os th, Z0 + (4 + 1)h. 
The s-th derivative of f(z) at any grid point z; can be expressed by 


(2) Ff (2) = Le Pe (z,)f (zx). 
Letting z = zo + Ph, it is possible to write 
(3) hf (23) = Die Mi (j)fz/ M(s). 


where f; = f(z), and M,“(j) and M(s) are complex integers, independent of 
f and h. The author tabulates the exact coefficients of f, for all grid points, in 
the derivatives of all orders, for polynomials of degree 2, 3, ---, 9. Methods 
of computing and checking are described in the paper. The work will add to 
the author’s reputation for supplying accurate and well planned tables of 
coefficients. 
GERTRUDE BLANCH 

NBSINA 


1H. E. Satzer, “Formulas for numerical integration of first and second order differ- 
ential equations in the complex plane, Jn. Math. Phys., v. 29, 1950, p. 207-216. 


1105[K ].—E. P. Kine, ‘‘The operating characteristic of the control chart 
for sample means,’ Annals Math. Stat., v. 23, 1952, p. 384-395. 


This paper extends the theory of the SHEWHART control chart by deriving 
expressions for the chances that the X chart or chart for sample means will 
show control for both of the cases of known and unknown standard devia- 
tion. The null-hypothesis is that samples are drawn from the usual normal 
process, N(u, o?), where uw and o are fixed but unknown. The alternative 
hypothesis considered is that shifts in the process mean from time to time 
can be represented by N (yu, @o?), i.e., the process mean itself is a random 
variable with this normal distribution. 

Suppose we consider m samples of m each which are selected at random 
from rational categories of a process and we plot the X and R (sample 
range) charts. Let Bo(k, 0, m, n) denote the chance that all m sample means 
will fall within +kon-! of the grand mean, where the standard deviation ¢ is 
known, and let 8(k, 0, m, n) represent the chance that the X chart will show 
control when ¢, in effect, is estimated from the m sample ranges on the chart 
for ranges, i.e., the process itself to date. Of course, it is usual practice to 
take k = 3 and the computed tables of the paper are based on this value. 
The Tables I-VI (p. 390-392) give values to 2D for the operating character- 
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istics or the chances, 8, that the X chart will show control for several practi- 
cal cases: Bo(3,0, m,n) for @ = 0(.5)3; nm = 2,5,10, m = 2,3,4; B(3,8, 
m,n) for @ = 0(.5)3; m = 2,5,10,m = 2, and m = 5,10, m = 3,4. 

For m > 4 bounds for 8o(3, 0, m, ) and 8(3, @, m, m) are given in Tables 
VII and VIII (p. 393, 394) to 2D for the cases: @ = 0(.5)3;” = 5,10; m = 
5,10; 6 = 0, .75, 1, 1.25, 1.5; 2 = 5,10; m = 15,20; 0 = 0, .25, .75, 1; 
nm = 5,10; m = 25. 

It is of interest to note that the problem treated in this paper is some- 
what related to that of testing for ‘‘outlying”’ observations, since the chance 
that the X chart will show control is also the chance that the largest and 
smallest sample means will both lie between control limits. 


F. E. GrRuBBs 


Ballistic Research Laboratories 
Aberdeen Proving Grounds, Maryland 


1106[K ].—W. H. Kruskat & W. A. WALLIs, ‘‘Use of ranks in one-criterion 
variance analysis,’’ Amer. Stat., Assn., Jn., v. 47, 1952, p. 583-621. 


Consider three samples, of sizes m1, m2, and 3 respectively, arranged to- 
gether in order of size. Assign scores to the N = nm; + m2 + nz individuals 
according to their ranks in the combined sample. That is, assign the score 
1 to the smallest of the N, the score 2 to the next smallest, etc. Let H = 
W-i° n{R; — 4(N + 1)F 

N z (NW? — 1)/12 
sample. 

Under the assumption that the three samples are random selections from 
a continuous population Table 6.1 (p. 614-617) of this paper gives, for m, 
N2, M3 < 5, exact probabilities (and three approximations to the probabili- 
ties) all to 3D, that H will equal or exceed certain selected values. The 
selected values are chosen so that the probabilities will be close to 10, 5, and 
1 per cent. 





, where R; is the mean rank score of the ith 


FRANK MASSEY 
University of Oregon 
Eugene, Oregon 


1107[K ].—L. E. Mosss, “‘A two-sample test,” Psychometrika, v. 17, 1952, 
p. 239-247. 


Samples of size m and n are drawn from populations A and B, respec- 
tively, and the combined samples are arranged in increasing order. To test 
the hypothesis that A = B against the alternative that B is more widely 
dispersed than A, the author proposes the statistic s,*, defined as one more 
than the difference in the ranks of the h-th largest and h-th smallest obser- 
vation in the first sample. He tables (p. 244, 245) the tail of the distribution 
of s,* to 3D for h = 3, 4,5, 8, m = 4h, and 5 suitably selected values of n 
in each case. 

J. L. Hopcgs, Jr. 
University of California 
Berkeley, California 
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1108[K ].—P. J. Riykoort, “‘A generalization of Wilcoxon's test,’’ Neder- 
landsche Akademie van Wetenschappen, Proceedings, v. 55, series A, 
1952, p. 394-404. 


Let x;; be the j-th observation in the i-th sample, where i = 1, ---, 
and j = 1, ---, ;. Let r;; be the rank of x,; in the combined sample of = 
>; observations. The author proposes testing the hypothesis that all of 
the samples come from distributions with the same mean value by the use 
of the statistic 


S= (ss; —n#)? 


where s; = }or;; and # = (m + 1)/2. 

The tables (p. 400-402) give to 3 and 4D (sometimes more) the exact 
cumulative distribution function of S: for k = 3 and each n; = 2, 3, or 4; 
k = 4 and each n; = 2; k = 5 and each n; = 2, for all values of S, and for 
k = 3 and each n; = 5; k = 4 and each n; = 3, for large values of S. 

Methods of approximating the distribution function of S by the chi 
square distribution and the analysis of variance distribution are presented. 
Using these methods a table (p. 402) is given showing the approximate five 
per cent points of S for all combinations of k = 3(1)10 and n; = 2(1)10. 

The statistic S is a linear function of the statistic H proposed by Krus- 
KAL & WALLIs.! These authors give tables of the five and one per cent points 
for H for k = 3 and all possible combinations of n; < 5. 

I. R. SAVAGE 
National Bureau of Standards 
Washington, D. C. 


1W. H. Krusxat & W. A. WaLLIs, “Use of ranks in one-criterion variance analysis,” 
Amer. Stat. Assn., Jn., v. 47, 1952, p. 583-621. (RMT 1106.] 


1109[(K ].—Co.tn Waite, “The use of ranks in a test of significance for 
comparing two treatments,”’ Biometrics, v. 8, 1952, p. 33-41. 


This paper presents tables for use of the WrLcoxon procedure! for com- 
paring two treatments when the numbers of individuals in the two groups 
are not necessarily equal. Let m; be the number of individuals in the group 
for which we compute the rank total T while mz is the number of individuals 
in the other group. Without loss of generality, m; can always be taken less 
than or equal to m2. The ranks allotted are 1, 2, ---, (a: + m2). The null 
hypothesis tested is that both treatments had the same effect ; that is, that 
T represents the sum of m,; ranks drawn at random from the finite universe 
1,2, ---, (m; + me). Let the integer JT, (notation of reviewer) have the 
property that Pr(JT < 7.) <a@ and Pr(T < T.+ 1) >a. Then also 
Pr[T > (mi + m2 +1) — Ta] <a and Pr[T > m,(m, + m+ 1) — Ta 
— 1] >a. Tables of T. (p. 37-39) are presented for a = 5%, 1%, .1% and 
all m;, m2 up to m, + m2 < 30. 

J. E. WatsH 
U. S. Naval Ordnance Test Station 
China Lake, California 


a. >; seaemae “Individual comparisons by ranking methods,”’ Biometrics, v. 1, 1945, 
p. ; 
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1110[L].—Le CENTRE NATIONAL pD’ETUDES DES TE£LE&COMMUNICATIONS, 
Tables des fonctions de Legendre associées. Editions de la Revue d'Optique, 
Paris, 1952. xxi + 292 p., 21 X 30.5 cm. 


Legendre functions whose degree is not an integer, arise in potential and 
wave problems relating to cones. Some computations of such functions have 
been reviewed previously (MTAC, v. 5, p. 152-153; v. 6, p. 98-99; see also 
RMT 1120) but no systematic tabulation seems to exist, apart from the case 
when the degree is an integer or half of an odd integer. Thus the present 
tables are a pioneering effort in a field which is becoming more and more 
important. 

The function tabulated in this volume is 


d™ P,, (cos 8) 


P,™ (cos 0) = (— sin 6)™ d(cos 6)" 


where m is a non-negative integer, @ is real, and 1 is real. 

The introductory material includes a Preface by J. CouLoMB; an Intro- 
duction by L. Rosin; an account of the method of computation, checks, 
accuracy by P. LE GALL; a diagram facilitating the use of the recurrence 
formulas; the level curves of P,° (cos 6) = const., ---, P,® (cos @) = const. 
in the region 0 < m < 10,0° < @ < 90° of the n, @ plane; and a detailed page 
index. 

The tables give values of P,,™ (cos @), to various degrees of accuracy, for 
m = 0 (1) 5,” = —.5 (.1) 10, 6° = 0° (1°) 90°. A second volume (in prepara- 
tion) will extend the tables to @ = 180°. 

A. E. 


1111[L].—P. C. CLemmow & Cara M. Munrorp, “‘A table of 4/ (42) et" 
e~t#™™* d) for complex values of p.”” R. Soc. of London, Phil. Trans., 


e 
v. 245A, 1952, p. 189-211. 


Tables of the error function for complex variable were long overdue. They 
are needed in many wave propagation, and in some other problems. The 
computation of this function was discussed in MTAC, v. 5, p. 67-70. 

The function chosen for tabulation in this memoir is 


G(p) = V (3m) ete? f e~ tir? dy. 


For large p with 0 < arg p < } = it is asymptotically represented by 


aeslt jh ng ee --| 
tr/ (29) p inp? (tmp)? 
and 

G(0) = 4rte-*#/4, 


The table gives 4D values of the real and the imaginary parts of G() for 
p = re” where r = 0 (.01) .8, 6° = 0° (1°) 45°. 


A. E. 
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1112[L].—H. M. Daccetrt, Jr., ““The Shedlovsky extrapolation function,”’ 
Amer. Chem. Soc., Jn., v. 73, 1951, p. 4977. 
4D tables of 
{32 + [1 + (42)? ]#}? 
for z = 0(.001).209. 
A. E. 


1113[L].—A. A. DoropniTfsyn, “Asimptoticheskie zakony raspredeleniia 
sobstvennykh znachenil dla nekotorykh osobykh vidov differenfsialnykh 
uravnenil vtorogo poriadka.’’ [Asymptotic laws of distribution of the 
characteristic values for certain special forms of differential equations of 
the second order. ] Uspekhi Matem. Nauk (N.S.) v.7,no. 6, 1952, p. 3-96. 
Two short tables (p. 95) give numerical values to 4-6 S of 


¢(s,a) = = (n + )-* 


Table I. s = 2(2)8, a = 3(1/12)4. 


Table II. s = 2(2)8, a = 3(.05)4. 
A. E, 


1114[L].—H. Gort er, “Zur laminaren Grenzschicht am schiebenden 

Zylinder. Teil I.’’ Arch. Math., v. 3, 1952, p. 216-231. 

The functions Fo, ---, F222 are solutions of the differential equations 
Fo” + fiFo’ = 0 
Fl’ + fiFe! — 2fi'F. = — 12 fsFo’ 

Fi’ + fiFd — 4fi' Fs = — 30 gsFo’ 
Foo + fiFo2 — 4f1'Fo2 = — 30 hsFo’ — 12 f3F 2’ + 8fs' Fe 
Fo’ + fiFe’ — 6f1'Fs = — 56 g7Fo’ 
Fo, + fiFoq — 6f1' Foe = — 56 hy Fo’ — 30 gsFo’ + 12 gs’Fo — 12 f3F,’ 
+ 16 fF 

Fase + f1Fo22 — 6f1' Foee = — 56 krFo’ — 30 hsF 2’ + 12 hs’ F2 — 12 fsF 22 

+ 16 fs’ Fee. 
All F’s vanish at 0, Fo = 1, and all the other F’s vanish at ~. The f;, g;, hi, 
k; are the functions appearing in the integration of the boundary layer 
equations according to BLastus! and Howartu.? 

The present paper gives 3D or 4D tables and diagrams of Fo(n), Fo’(n), 
-++, Feoo(m), Fae2(n) for 7 = 0(.1)5.4 (except for Fos, ---, F222 when » < 
5.2). The values have been computed by numerical integration of the differ- 
ential equations. The necessary values of f;, ---,%; were taken from the 
tables by ULricu,' and subtabulated where necessary. The values of Fy and 
F,’ were compared with those given by Cooke‘ and ScHLICHTING;' dis- 
crepancies in the last decimal were noted. 


A. E. 


1H. Brastus, “Grenzschichten in Fliissigkeiten mit kleiner Reibung,” Zschr. f. Math. 
u. Phys., v. 56, 1907, p. 1-37. 

2L. HowartH, Steady flow in the Boundary Layer near the Surface of a Stream. ARC 
Report No. 1632, 1934. 
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3 A. Utricu, “Die ebene laminare - an einem Zylinder,’’ Arch. d. Math., 
v. 2, 1949, oP 37-41; MTAC, v. 4, 1950, p. 96-97. 
OOKE, “The boundary layer of a class of infinite yawed cylinders,’’ Cambridge 
Phil, , Proc., v. 46, 1950, p. 645-648. 
5H. SCHLICHTING, Grenzschicht-Theorie. G. Braun, Karlsruhe, 1951. 


1115[L].—J. L. Lusxin & Y. L. LuKE, “Frequencies of longitudinal vibra- 
tion for a slender rod of variable section,” Jn. Appl. Mech., v. 20, 1953, 
p. 173-177. 


The frequency equation is 


[Ji(a) ¥1(8) — Yila)Ji(8)] sin y + (€/le|)[ Ji (a) Yo(8) 
— ¥i(a)Jo(8)]cos y = 0 


where a = (1+ €)8, B = vaLe/le|, y = vali and », v2---are proportional 
to the frequencies. The dimensionless parameters ¢ and A = L,/(L; + Lz) 
are used, and table 1 gives 5 D values of »v,(Z: + Lz) for m = 1 (1) 5, 
e = — 1 (1/3) 1, X = O (.125) 1. “Seven or eight decimals are carried in 
the computations, sufficient to insure that only occasional entries, if any, 
are in error by as much as one unit in the fifth place.”’ 

Table 2 gives some auxiliary quantities, useful in interpolating in the 
e-direction. 

A. E. 


1116[L].—D. MANTERFIELD, J. D. CRESSWELL, & H. HERNE, “The quick- 
immersion thermo-couple for liquid steel,’’ Iron and Steel Institute, Jn., 
v, 172, 1952, p. 387-402. 


Table I, p. 396, gives 4D values of the first six roots x of the equation 
Jo(kx) Yi(x) — Yo(kx)Ji(x) = 0 
for k = (1.1)? 1.06(.02) 1.1(.05) 1.3(.1) 1.5, 2(1), 5. 
Table II gives the corresponding 4D values of 
2F1 (x) Jo(Rx) 
x{[Jo(kx) P — [Ji (x) P} 
for k = (1.1)#, 1.15(.05) 1.3(.1) 1.5, 2(1) 5. 


The tables were computed by H. Carsten and N. McKeErrow;; they are 
referred to in FMR Index, section 17.812, p. 268. 





A. E. 


1117[L ].—NBS, Applied Mathematics Series No. 13, Tables for the Analysis 
of Beta Spectra. U. S. Government Printing Office, Washington, D. C., 
1952. iii + 61 p., 20 X 26.5 cm. $0.35. 

These tables were designed to assist in the theoretical analysis of beta- 


ray spectra. 
The principal tables give values of the “‘Fermi function” 


F(Z, n) = 5 et |P(1 + S + 46)/? 
where Z is the atomic number, 7 the momentum of the electron, 
> (1 + 7’)}, 8 = n/€,¥ = 2/137, S = (1 - 7)! _ 1,6 = 7/B. 
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Other abbreviations used in this pamphlet are 
T = 510. 91(e — 1), Bp = 1704 .3n. 


A brief introduction gives definitions and indicates the method of com- 
putations, and is followed by sections on beta spectra and their analysis, 
and a bibliography. 

There is a set of 6 auxiliary tables. 

Table 1 (p. 15-16). Values of 9, ¢, 8, T[4S or 5S) for Bp = 0(100)20000 
Gauss cm. 

Table 2 (p. 17). 5S values of ¢, 7, 8, Bp for T = 0 (10) 60 (20) 200 (50) 
1000(100) 5000 (200) 10000 kev. 

Table 3 (p. 18). Values of the “nonrelativistic” electrostatic effect factor. 

Table 4 (p. 18). Illustrating the importance of the relativistic effect for 
Z = 0, 10 (20) 70, 100 and » = 0, .2, .6, 1, 2, 4, 7. 

Table 5 (p. 18). Values of the ratio 





ri+S+i)|.. as 
for the same values of Z, 7 as in Table 4 except that 7 = 0 is omitted. 


1 a vit+# 
5.109dVi+_ 7 
and of ae , to 4D (two separate tables) for Z = 0, 10, 20 (20) 100, » = 
0 (.2) 1 (1) 6. 

Note that the ‘‘Fermi-function”’ f(Z, 7) is not identical with the ‘“Fermi- 
Dirac function” F, (7). 

The principal table (p. 21-61) gives values of f(Z, ») for Z = 1 (1) 100, 
n = 0 (.05) 1 (.1) 7. Each Z has a column to itself going over two adjacent 
pages. The column heading gives values of Z, y, S, ¢(Z), and each row (for 
a given n, T, ¢, Bp all of which are stated) gives the values of f(Z, 9) indi- 
cated as B- and 8+ for the upper and lower signs in the exponential respect- 
ively. At the foot of each column constants A,B,C are listed, and computation 
from the approximate formula 


BC 
(2,9) = 8(4 +245) 


is recommended for 7 = 0. 
On p. 10, and again on p. 21, ¢(Z) is defined as 


(4rmcR/h)*S(T (3)/T (3 + 2S) P(1 + S/2), 


and this definition was used in the auxiliary tables. In a mimeographed cor- 
rection sheet (which should accompany each copy) it is pointed out that the 
values listed as g(Z) on p. 22-61 of the principal tables are actually values 


of 
(EX 2 





Table 6 (p. 18, 19). Values of f(Z,) to2 or 3S 





3+ 53) (2). 


This sheet corrects the value for Z = 100, and lists 5S values of ¢(Z) for 
Z = 1 (1) 100. 


A. E. 








182 RECENT MATHEMATICAL TABLES 


1118[L].—NBS, Applied Mathematics Series No. 25, Tables of Bessel 
functions Yo(x), Yi(x), Ko(x), Ki(x), O0< x <1. U. S. Government 
Printing Office, Washington, D. C., 1952. ix + 60 p., 20 X 26.5 cm. 
$0.40. 


This is a reprint of Applied Mathematics — No. 1 and was reviewed 
in MTAC, v. 3, p. 187-188. 
A. E. 


1119[L].—NBS, Applied Mathematics Series, No. 28, Tables of Bessel- 
Clifford functions of orders zero and one. U. S. Government Printing 
Office, Washington, D. C., 1953. ix + 72 p., 20 X 26.5 cm. $0.45. 


“Although the Bessel-Clifford functions are obtainable from existing 
tables of Bessel functions, it was felt that they warranted tabulation because 
it is generally necessary to enter the existing tables with an irrational argu- 
ment. Furthermore, the Bessel-Clifford functions arise as solutions of a class 
of differential equations occurring in various branches of applied physics, 
and they are therefore of importance in themselves.” “‘The present volume 
carries the tabulation of the functions of orders zero and one up to a point 
where the asymptotic expansions can be used conveniently.’’ (From the 
introduction.) 

Table I. Jo(2yx) for x = 0 (.02) 1.5 (.05) 3 (.1) 13 (.2) 45 (.5) 115 (1) 
410, 8D. J,(2-Vx)/ Vx for x = 0 (.02) 1.5 (.05) 3 (.1) 13 (.2) 45 (.5) 115 (1) 
125, 8D, x = 125 (1) 410, 9D. 

Table II. Y,(2x)(x)-", same values of m and x as in Table I. 

Table III. I9(2-¥x) for x = 0 (.02) 1, 8D, x = 1 (.02) 1.5 (.05) 6.2, 7D. 
I, (2¥x)/vx for x = 0 (.02) 1.5 (.05) 6.2, 7D. 

Table IV. e~*¥? I9(2Vx) for x = 6.2 (.1) 13 (.2) 36 (.5) 115, 8D, x = 115 
(1) 160 (5) 410, 9D. e°¥* I,(2¥x)/ Vx for x = 6.2 (.1) 13 (.2) 36 (.5) 65, 
8D, x = 65 (.5) 115 (1) 160 (5) 410, 9D. 

Table V. Ko(2Vx) for x = 0 (.02) 1.5 (.05) 2.5, 8D, x = 2.5 (.05) 6.,2 
9D. K,(2Vx)/ Vx for x = 0 (.02) .04, 6D, x = .06 (.02) .28, 7D, x = .30 
(.02) 1.5 (.05) 2.5, 8D, x = 2.5 (.05) 6.2, 9D. 

Table VI. e?¥* Ko(2-x) for x = 6.2 (.1) 13 (.2) 36 (.5) 115 (1) 160 (5) 
410, 8D. e”* K,(2yx)/Vx for x = 6.2 (.1) 13 (.2) 36 (.5) 40, 8D, x = 40 
(.5) 115 (1) 160 (5) 410, 9D. 

Except in the neighborhood of singularities, second central differences, 
sometimes modified, are given in all six tables. 

Auxiliary tables. Interpolation coefficients for Everett’s formula. 

The Introduction (by M. ABRAMOWITZ) gives the mathematical proper- 
ties of Bessel-Clifford functions, describes (mathematical) applications, in- 
terpolation, the method of computation, and accuracy, and lists references. 

A preliminary tabulation of these functions was reported in MTAC, 
v. 3, p. 107. 

A. E. 
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1120[L ].—K. M. Srecet, J. W. Crispin, R. E. Ktemman, & H. E. Hunter, 
“The zeros of P,, (xo) of non-integral degree.” Jn. Math. Phys., v. 31, 
1952, p. 170-179. 


Let P{?(x) be the associated Legendre function of order unity and non- 
| integral degree. This paper is concerned with determination of m; and 


1 
f [PS (x) ? dx such that P{?(x) = 0. The evaluation of m; has previously 
z0 


been considered (MTAC, v. 5, p. 152-153; v. 6, p. 98-99). The idea of the 
paper is to expand P&?(x), n; = n+ Z,, in a Taylor series and calculate 
Z, by inversion. The formula for Z, depends on Legendre and associated 
Legendre functions and a certain sum. The latter is finite if m is an integer or 
an odd multiple of 1/2, and for these choices of computation is facilitated 
as many tables exist for the former elements. The series expansion is also 
used to evaluate the above intergral. 

The theory is illustrated for x» = cos 165°. In Table 3, values of P.(xo) 
and P§? (x9)/V1 — x0? are tabulated to 10S for m = 1 (1) 20. Values are also 
given for m = — 0.5(1.0)20.5; the number of significant figures varies from 
5 to 11. These were calculated using series expansions and the authors state 
that existing tables were used as checks. The entries for 2 an odd multiple of 
1/2 are mostly new. For m an integer, many of the entries are already avail- 
able (cf. FMR Index, for example). 


In Table 4, the first 19 values of m; and f [PS (x) dx are given. 
=z0 


Three terms of the Taylor series are used. In each case computations are 
presented for both m an integer and an odd multiple of 1/2. No estimates of 
the error are given, but a procedure is defined to indicate preferred value 
according to the choice of ”. The authors are interested in a speedy means to 
obtain entries within engineering accuracy, and the procedure outlined 
seems to fulfill this need. 

YupELL L. LUKE 
Midwest Research Institute 
Kansas City, Missouri 


1121[L].—Frieprich TOLKE, Praktische Funktionenlehre. Erster Band. 
Elementare und elementare transzendente Funktionen. Springer-Verlag. 
Berlin, Géttingen, Heidelberg. 1950. xii + 440 p., 20 27.5 cm. 


The first edition of this volume appeared in 1943 and was not accessible 
to the reviewer. The present second edition is said to be considerably en- 
larged. Further volumes, on theta- and elliptic functions, hypergeometric, 
Bessel, and Legendre functions are in preparation. The whole work will 
present the more important functions in a manner suitable for their applica- 
tion in the engineering sciences. 

The present volume contains a number of numerical tables. They are 
designed for use by engineers rather than professional computers. There are 
instructions for using these tables, but no sources or indications of the method 
of their computation, or of their accuracy. 

Table 1 (p. 168-207), 2x, In 2xx, e**, e***, sinh 24x, cosh 2x, 
tanh 2x, coth 2xx, amp 2x, sin 24x, cos 24x, tan 24x, cot 24x, sin* 2rx = 
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sin (24x — 39), cos* 2x, tan* 2x, cot* 2x, to 4 or 5D, mostly with first 
differences, for x = 0 (.001) 1. 

Table 2 (p. 208-247). 4ax, e!**, e~#**, sin 44x, cos 42x to 4 or 55S, for 
x = 0 (.01) 40. 

Table 3 (p. 248-257). Ei(x), Ei(—x), Shi(x), Chi(x), Si(x), Ci(x) to 
3-4S, with differences, for x = 0 (.01) 5. Here 


Bi(e«) = In |x| + a eae = Chi(x) + Shi(x) 
Silx) = — ¢Shi(éx), Ci(x) = Chi(éx). 


Auxiliary tables on a folding inset to p. 258. 4x, e!**, e~#**, sin 42x, 
cos $ax, 3-5S, x = .001 (.001) .01. 2ax, e**, e-***, x = 0 (1) 10. e”, e-, 
x = 1 (1) 50. 
marx mx kax ax 
td’ 3° 3k 5D, m = 1 (1) 8, k = 1 (1) 4. 


m!, 1/m!, m = 2 (1) 10. m!/(m —n)!, n =1(1)m, m=1(1)10. (), 
m = 1 (1), = 1 (1) 15. Some useful numerical data. 
Functions which occur in diffusion problems. 


Fy(x, y) = 4 ez cos 2kry, | 

5D, x = 0(.01).05(.05).25 (.25)1.5(.5)2.5, 4, y = 0(.01).5(p. 268-9). 
F,(x, y) = x (kn)-! e-"* sin 2key, 

4D, x = 0(.01).05 (.05).5(.1)1(.2)2, 2.4, 2.8, y = 0(.01).5(p. 271-2). 
F;(x, y) = z + x (—1)*(2e°k?)"[1 — e-***] cos 2kry, 

4D, x = 0(.01).05(.05).5(.1)1(.2)2, y = 0(.01).5(p. 275-6). 

F,(x, y) = x (4x*k)1[1 — e-*=] sin 2key, 


4D, x = 0(.0025).0125 (.0125).125 (.025).25 (.05).5(.1)1(.25)2.5, 
y = 0(.01).5 (p. 282-4). 
Legendre polynomials P(x) m = 0 (1) 10, 4D, their derivatives P,,’ (x), 
n = 3 (1) 6, 4D, and their integrals 
Pix (x) - f ~— f P,,(x) (dx)*, 5 = 0(1)6, k= 1, 2, 5D, 
0 0 


for 
x = 0(.001)1 (p. 372-440). 


The volume also contains tables of integral formulas for indefinite 
integrals (p. 69-156). 


A. E. 
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41 


101 
103 
127 


25 
49 
121 





1877, p. 86-88. 

08741 

26380 36360 

05135 

38410 86376 34932 
93399 

63766 71392 

76679 

80965 


The above list of errata is complete. 


E.M.P.T. 
Tulle (Corréze) 


France 


1940, p. 123. 


for 
399 
108 
830 
756 
921 
445 
829 
183 
051 
233 
702 
833 
163 
224 
658 
005 
127 
024 
767 
742 


087 
984 
419 
282 
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In this issue reference has been to errata in RMT 1097. 


38409 


26379 


74238 


63763 


08739 


82105 
05185 
64260 
95958 


48908 
76675 
81074 


227.—S. M. Dracu, “‘Cube roots of primes to 31 places,"’ Mess. Math., v. 7, 


read 
228 
110 
860 
760 
915 
726 
830 
597 
119 
790 
893 
789 
148 
226 
768 
948 
269 
023 
765 
743 


308 
923 
408 
283 


R. Li£NARD 


228.—(a) F. Empe, Tafeln Elementarer Funktionen, Leipzig and Berlin, 


(b) M. Bott, Tables Numériques Universelles, Paris, 1947, p. 472. 


Both these works contain a table of Langevin’s function 


L(x) = cothx — x“. 
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Comparison of these with a recently computed unpublished table (see 
UMT162) reveals the following errata. 


x EMDE BoLL L(x) 
0.20 .06649 .067 .066490 
.30 .09941 .100 .099405 + 
.36 11897 — 118976 
48 15759 — -157595+ 
.60 .19536 .196 .195359 
72 .23210 — .232095 — 
1.36 4058 — .405746 
1.72 4848 — 484858 
1.86 5121 — .512037 
1.90 .5195 .520 .519450 — 
2.08 .5510 — .550941 
2.50 .6135 .614 .613567 
2.60 .6265 .624 .626479 
3.65 7273 — .727379 
3.70 .7309 re | .730953 
4.00 -7507 .749 .750671 
4.20 .7623 .762 .762355 — 
4.50 .7780 177 .778025 — 
8.20 .8781 — .878049 


ANDREW YOUNG 
Univ. of Liverpool 
Liverpool, England 


229.—R. A. FISHER & F. YATES, Statistical Tables for Biological, A gricul- 
tural and Medical Research, 3rd ed., New York, 1948. 


Table VI gives significance levels for the correlation coefficient. Errors in 
the table of the variance ratio previously reported (MTAC, v. 6, p. 35-38) 
have affected Table VI for » = 0.001 as follows: 


n For Read 
3 0.99116 0.99114 
5 0.95074 0.95088 

45 0.4648 0.4647 

70 0.3799 0.3798 


In addition, an error has been introduced at m = 6, where 0.92493 should be 
0.92490. There are also a few ‘rounding errors.’’ Two of these are worth 
recording because they constitute disagreements with the corresponding 
values of z. For nm = 13, r = 0.5139 should be 0.5140, and for nm = 45, r = 
0.2875 should be 0.2876, both for p = 0.05. 

H. W. Norton 


University of Illinois 
Urbana, Illinois 
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ee 230.—K. Hayasul, Tafeln der Besselshen, Theta-, Kugel- und anderer Funk- 
tionen, Berlin, 1930. 
Table X of 12D values of square roots contains the following errata. 
The starred entries have been noted by Fukuoka. There are 27 other 
errors in the 11th and 12th places. 
For Read 
0,217 0,46583 15146 06 0,46583 25879 54 
0,289 0,53758 72032 29 0,53758 72022 29* 
0.498 0,70554 94312 95 0,70569 11505 75* 
0,543 0,73688 53302 92 0,73688 53370 78* 
0.699 0,83606 21926 13 0,83606 21986 43 
0,787 0,88713 02040 49 0,88713 02046 49 
0,800 0,89442 70910 00 0,89442 71910 00 
0,845 0,91923 88101 03 0,91923 88155 43 
0,917 0,95760 11173 76 0,95760 11695 90 
2,53 1,59057 73720 59 1,59059 73720 59 
2,90 1,70293 86395 93 1,70293 86365 93 
3,17 1,78044 93812 80 1,78044 93814 76 
6,40 2,52982 20281 35 2,52982 21281 35 
7,36 2,71293 18932 50 2,71293 19932 50 
8,55 2,92403 83031 93 2,92403 83034 43 
8,71 2,94127 09126 75 2,95127 09126 75 
8,72 2,94296 46120 47 2,95296 46120 47 
8.73 2,94465 73405 39 2,95465 73405 39 
8.74 2,94634 90998 19 2,95634 90998 19 
8,75 2,94803 98915 50 2,95803 98915 50 
8.80 2,96647 93848 38 2,96647 93948 38* 
8.98 2,99666 28127 54 2,99666 48127 54* 
ul- 9,31 3,05122 92504 78 3,05122 92604 87 
9,51 3,08382 87889 22 3,08382 87890 22 
™ 14,1 3,75499 66666 56 3,75499 66711 04 
38) 14,3 3,78153 40798 26 3,78153 40802 38 
FRIEDEMANN SINGER 
Robert-Blum-Strasse 11 
Halle (Saale) Germany 
Lupwic STAMMLER 
Reilstrasse 129 
Halle (Saale) Germany 
- UNPUBLISHED MATHEMATICAL TABLES 
rth 159[A].—NATIONAL PuysicaL LABORATORY (Great Britain), Tables of Bi- 
ing nomial Coefficients. 20 quarto pages. Deposited with the Roya Society 
= (no. 15). 
The Binomial coefficients C, = are given to six decimal places for 
x = 0(.001)1; m = 2(1)8. These are coefficients in the Newton-Gregory 
formula for interpolation with forward differences. 
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160[D,L].—InstitutT pE MatuémaTigues AppLiguf&es, Ecole Polytech- 
nique de |’Université de Lausanne (Switzerland). Manuscript in posses- 
sion of the author. 


Table of the function F(x) = sin x/x and of its derivatives d"F(x)/dx" 
for x = 0(.01)4, m = 0(2)16, 10D. : 
This table has been obtained by subtabulation from another table cal- . 


culated on the base of series for x = 0(.1)4.3. A National accounting ma- t 

chine (class 3000) was used for the subtabulation. The error in the original 

tables does not exceed .55 units of the last order. m 
Cu. BLANC I 

161[D ].—SuBMarRInE SIGNAL DIVISION OF THE RAYTHEON MANUFACTURING 1 


Company, Transducer Department, Boston, Sines and Cosines of the 
Decimal Circle. 10 lithographed octavo pages of tables. Deposited with 
the Roya. Society (no. 2). 


This table gives sina and cosa for a/2x = 0(.001)1 to 5D when the 
entry is less than 1/6 and to 4D otherwise. 


162[E ].—A. Younc & T. Murpny, “Tables of the Langevin Function and 
its inverse,” 6 leaves. Deposited in UMT File, also in the depository of 
the Roya Society (no. 21). 


The function L(x) = coth x — x is tabulated for x = 0(.01)7.50 to 
6D. 

A comparison of this table with previously existing tables reveals errata 
in the previous tables which will be found in MTE 228. 


163[G ].—E. M. IBraui, Tables for the plethysm of S-functions. Deposited 
with the Roya. Socrety (no. 1). 


This comprises about two dozen large sheets of manuscript. The purpose 
of these specialized algebraic tables is described by the author in Quar. Jn. 
Math. (Oxford Series), Ser. 2, v. 3, 1952, p. 50-55. They relate to a special 
type of multiplication connected with symmetric functions and extend to 
partitions of a total degree of 18. 


164[1,L ].—C. W. Jonss, J. C. P. Miter, J. F. C. Conn, & R. C. PANK- 
HURST, Tables of Chebyshev Polynomials. 2 typed double-foolscap sheets. 
Deposited with the Roya Society (no. 7). 


These tables give exact values of C,(x) = 2 cos(mcos—'4x) for x = 
0(.02)2 and m = 1(1)10. For m = 8, 9, 10 they extend and complete the 
curtailed values given in the table described in RMT 381, MTAC, v. 2, p. 
262. [See also RMT 1103. ] 


165[L].—ApDMIRALTY RESEARCH LABORATORY, Solution of the Equation 
(y’’)? = yy’ and two other Equations. 3 foolscap manuscript pages of tables 
and 4 pages of description. Deposited with the Roya Socrety (no. 9). 


The solution of the equation (y’’)? = yy’ for which y = 0, y’ = 1 when 
x = 0, is tabulated to 6 decimals for x = 0(.05).5(.1)6 and facilities are 
provided for interpolation. 
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The integral 


r) 
22-38 exp( — p’) f eIo(2pe) exp( — e)de 


is tabulated to 4 decimals for the range 8 = 0(.25)4, p = 0(.25)5. Interpola- 
tion is possible in this table, but no differences are provided. 

The root x of the equation u sin x — cos x + e~“* = 0 that lies between 
x and 2z is tabulated to 4 figures for u = .1(.01).3(.02)2 and u# = 0(.02).5. 
Interpolation is linear. 


166[L].—M. S. Corrincton, Tables of Fresnel integrals, modified Fresnel 
integrals, the probability integral, and Dawson’s integral. Radio Corpora- 
tion of America, R.C.A. Victor Division. 25 quarto pages. Deposited with 
the Roya Society (no. 4). 


These tables give values for x = xu? = 0(.001).02(.01)2 of the functions 


C(u) = 3 f * T4(Oat = f ” cos (}af*)dt 
S(u) =} f * Jy(t)dt = f " sin (3af*)dt 
Ch(u) =} f “(dt = f " cosh (}xf)dt 


Sh(u) = 4 f " Ty(Odt = f " sinh (}xf?)dt 


v2 “ z= et 2 zi . . 
H (x!) = =f Ki; (t)dt = f et = =f, e~" dt 


and 


v2 —z = et 2 zt - 
D(x!) = — , Kyat = f nat = =f e* dt 


Two versions of the tables are given, one to 5D and another to 8D, with 
an error up to 2 final units. 


167[L ].—C. Mack & M. Cast ez, Tables of 


f Io(x)dx and f Ko(x)dx. 
0 a 
Deposited with the Roya Society (no. 6). 


The integrals have been tabulated to 9D for the range of argument 
a = 0(.02)2(.1)4. A brief description of the method of computation is given 
and also of the extent to which the tables are interpolable. 
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168[L ].--NATIONAL PuysicaL LaBoratory, Tables of the complex Jacobian 
zeta function. 9 foolscap pages. Deposited with the RoyaL Society (no. 
14). 
The Jacobian Zeta function Z, of modulus k = sina can be found from 
the tabulated function f, by the relation 


Z,(Ky1 + 1K'$,k) = fi (V1,6,@) -+- tf2(W1,6,a) > sirg/K 
f2(W1,0,a) _ fi(l = V1, 1 — ¢, 4a = a) 
where K, K’ are the complete elliptic integrals of modulus k. f; to 3 significant 


figures is given for ¥: = 0(.1)1; @ = 0(.1)1; a = 5°(5°)85°. No provision is 
made for interpolation. 


169[L ].—NATIONAL PuysicaL LaBoraTory, Integrals of Bessel Functions. 
2 quarto pages. Deposited with the Roya. Society (no. 17). 


10D values of f Jo(t)dt and f Yo(t)dt for x = 0(.5)50. No provision 
0 0 


is made for interpolation. 


170[L ].—NATIONAL PuysicaL LaBoratory, Table of 


Qe 
f J (2k sin 30) cos? $6 dé. 
0 
1 quarto page. Deposited with the Roya Society (no. 18). 


4D values are given for k = 0(.1)10. The table is interpolable using 
second differences, but no differences are given. 
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Edited by the Staff of the Machine Development Laboratory of the National Bureau 
of Standards. Correspondence regarding the Section should be directed to Dr. E. W. 
Cannon, 415 South Building, National Bureau of Standards, Washington 25, D. C. 


TECHNICAL DEVELOPMENTS 
A SPECIAL PURPOSE DIGITAL COMPUTER 


1. Design considerations. The design considerations of a special purpose 
computer for the solution of a large number of simultaneous linear, algebraic 
equations depend not only on the number of equations with which the com- 
puter must deal but also upon the properties of the matrix of the equations, 
the time to be allowed for computation and the required accuracy of the 
solution. In this particular case, it was estimated that up to 1200 equations 
might be expected and the arbitrary time of one day was allowed for com- 
putation, after the problem had been set up on the computer. An accuracy 
of one part in 100 was demanded in the solutions. 

A set of simultaneous equations can be represented in the matrix notation 
as 


(1) AX =C 
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where A is the matrix of coefficients, a,;;, X is the column vector of unknowns 
x,;,and C is the column vector of constant terms on the right hand side of the 
equations, c;. The first equation in the set of equations would look like this. 


(2) 1X1 + Gy2x%X%2 +°++ + inXn = C1. 


There are many methods of solving such equations ' but the most suit- 
able, from the standpoint of machine computation, is that of Gauss-Seidel.” 
This procedure is an iterative procedure, each successive iteration yielding 
a better estimate of the desired vector X. Thus at the end of the first iteration 
one has X®, at the end of the second, X, etc. One is assured that X‘” ap- 
proaches X as n becomes very large if the equations are normal equations. 

The computation associated with the Gauss-Seidel procedure is as 
follows: 


(3) xP = (¢; — aaxi® — ajaxe® +++ — a, eatin 
— a x) «mee *dinXn?) /ais 


where x” is the i** unknown being calculated in the p* iteration. Equation 
(3) simply states that to compute x, one takes c;, subtracts the products 
indicated, and divides the whole by a;;. The products in eq. (3) consist of 
two kinds, first products involving x; and x,-». That is, as an unknown is 
calculated in the p* iteration, this new value of the unknown is used in suc- 
cessive calculations in that iteration. 

Two different methods of performing the computation indicated in eq. 
(3) suggest themselves. The difference lies in the method of storage. One may 
feed the coefficients, a;;, into the arithmetic section of the machine according 
to rows and the unknowns, x,;), in synchronism with the coefficients and 
perform the desired operation. The diagonal terms, a;;, must be left at the 
end of each row as they are the divisor of the accumulated sum. This method 
requires that each coefficient be included in order that proper synchronism 
be maintained, even though the coefficient be zero. The alternate to this 
method consists of feeding only non-vanishing matrix coefficients into the 
arithmetic section, as before, and storing the unknowns in such a fashion 
that they can be called from storage in an arbitrary order. This allows con- 
siderable saving in storing matrix coefficients, if the percentage of non-van- 
ishing terms is small, at the expense of some complexity in the storage of the 
unknowns. In the problems we expected to deal with, only 1% of the co- 
efficients were of the non-vanishing variety and the second method outlined 
above was selected. Fig. 1 is a block diagram of the computer using this 
method. 

2. Storage. For a given problem, the order of the matrix coefficients 
entering the computation is fixed and furthermore is periodic, with period of 
one iteration. Magnetic tape was selected as the medium for storing these 
terms, a continuous loop of tape being used. The matrix coefficients are re- 
corded serially along one channel of the multi-channel tape, along with an 
order that specifies the type of term, i.e., ¢;, a4; or a;;. In an adjacent channel 
of the tape is recorded the column address of the matrix coefficient. It is this 
address that is used to specify the location of the corresponding unknown, 
x;. The coefficient-order, called a word, is stored in 25 binary digits. A front 
view of the computer is given in the frontispiece. This shows the tape trans- 
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port mechanism with a continuous loop of tape. The tape is run at 15 in./sec. 
and the coefficients are read off at 20 words per second. 

Since the unknowns do not appear in the computation in a systematic 
fashion, it is necessary to provide random access to this storage, with the 
access time commensurate with the speed at which matrix coefficients are 
produced. A magnetic drum was selected for this storage. The drum revolves 
at approximately 59 rps and thus maximum access time to any word on the 
drum is roughly 17 milliseconds. The capacity of the drum storage is dictated 
by the maximum number of equations to be handled; in this case one must 
allow for storage of 1200 words. Fig. 2 is a photograph of the magnetic drum. 

Words on the drum are stored serially, 64 words of 25 binary digits in 
each of the 19 tracks around the periphery of the drum. A ‘‘return to zero” 
system of pulse recording is used and pulses are recorded at a density of 
about 60 pulses per inch. The desired track on the drum is obtained by 
switching a matrix of gates. This switching is done by inserting the track 
address into an address register, AR. The desired angular position on the 
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Fic. 1. Block Diagram of Computer 


drum is obtained by counting the number of words that have passed under 
a head, and noting coincidence between this count and the angular position 
as specified by the AR. When coincidence occurs, the operation (read or 
write) may proceed. This BOW (Beginning of Word) counter is set to zero 
once every revolution by an origin pulse. 

3. Arithmetic section. All computation is performed in the arithmetic 
section of the computer. The following operations and the time required for 
each are indicated below. 


Addition 200 microseconds 
Subtraction 200 microseconds 
Multiplication 8 milliseconds 
Division 8 milliseconds 


The operation times listed above obviously do not include 17 milliseconds 
(maximum) required to secure one of the operands from the drum storage. 
All computation is performed in parallel,? that is the four operations 
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above are performed as parallel transfers from register to register. Carries 
proceed by means of delay circuits between digits. Thus for an m digit 
register, nd seconds must be allowed for carries to occur, where d is the delay 
time between digits. No time advantage over serial operation is obtained 
but it is possible to modify this method of operation to include a simultane- 
ous carry* when faster storage becomes practical. 

There are three registers in the arithmetic section, a multiplier register, 
MR, multiplicand register, MD, and a product accumulator, PA. The fol- 
lowing transfers are possible: + MR to PA, + MD to PA. Of the 25 binary 
digits in the MR, MD registers, 20 digits are occupied by numbers and 5 
digits by order. The binary point may be located either to the right of the 
least significant digit, or to the left of the most significant digit. The PA 





Fre. 2. Magnetic Drum Storage Unit 


register contains 40 binary digits and the binary point is located at the center 
of the register. Each register may be used as an accumulator or as a shifting 
register. 

In normal computation the matrix coefficients are read from the mag- 
netic tape into the MR register, the address of the coefficient into AR, the 
corresponding unknown is read into the MD and the indicated operation is 
performed. In the case of division, the quotient x” is formed in the MD and 
is compared, digit by digit, with the number previously in the MD, i.e. 
x?-), Thus too great a difference in the two values indicates an error in 
computation (and the new value of the unknown is rejected) and very minor 
differences for all the unknowns may serve as a criterion for convergence of 
the solutions. 

The frontispiece is a photograph of the complete computer. The rack on 
the left contains the power supply. The two racks adjacent to the power 
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supply carry the three registers of the arithmetic section and associated con- 
trols. The next rack from the left carries the address register, AR, and drum 
controls. The tape rack on the right has already been discussed. 

4. Logical design. The computer is a purely binary machine. Numbers are 
represented in the binary notation in order to take advantage of the sim- 
plicity of arithmetic operations with this notation. Negative numbers are 
represented as the “‘one’s’”’ complement of the absolute value of the number. 
That is, if the number 10 is written as 0000001010, the number —10 is 
written as 1111110101. In other words, the complement of the absolute value 
of the number has been taken with respect to the modulus of the register 
minus one, M — 1. With this representation of negative numbers, the prod- 
uct of two numbers may be formed without taking cognizance of the signs 
of the numbers as such. Let —A be represented as C(A), the complement 
of A. Then 

C(A)-B = C(AB) 
and 


C(A)-C(B) = AB 


if the partial products are accumulated modulo M — 1. The modulus of an 
accumulator may be reduced by 1 by the addition of ‘“‘end around carry.” 
It consists of feeding a carry from the last digit of the register back to the 
first digit. It is easy to show that the modulus of a digit register, modulus 
2”, is reduced to a modulus of 2” — 1, with the addition of end around carry. 

5. Logical elements. The logical elements or components of the computer 
were designed on the basis of obtaining reliable operation over wide varia- 
tions of supply voltages and tube characteristics. The components were built 














Fic. 3. Typical Plug-in Unit 
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up on plug-in chassis, as shown in Fig. 3, for ease of maintenance of the com- 
puter. The components include flip-flops, delay circuits, gates, and cathode 
followers. Each element requires a 15-20 volt positive pulse input and its 
output may be of either polarity. The elements are insensitive to negative 
pulse inputs. Around 400 of these elements are used in the computer. 

6. Operation. The computer has been in operation about 3 months and 
two problems have been set up and solved. One consisted of 73 equations 
with a computation time of 3 minutes per iteration. The other consisted of 
793 equations and 20 minutes of computing time were required per iteration. 
In both cases, 15 iterations were sufficient to yield solutions of the desired 
precision. 

J. P. WALKER, Jr. 
Haller, Raymond and Brown, Inc. 
State College, Pa. 


1 FRAZER, DUNCAN & CoLiar, Elementary Matrices. Cambridge, 1946. 
2 WHITTAKER & Rosinson, The Calculus of Observations. London, 1924. 
. 3 ENGINEERING RESEARCH AssOcIATES, High Speed Computing Devices. McGraw Hill, 
1950. 
4 INSTITUTE FOR ADVANCED Stupy, Preliminary Discussion of the Logical Design of an 
Electronic Computing Instrument. Princeton, 1947. 


BIBLIOGRAPHY -OF CODING PROCEDURE 


The material described below is among that which has been added to 
the collection at the National Bureau of Standards Computation Labora- 
tory. An equivalent collection is available at the National Bureau of Stand- 
ards Institute for Numerical Analysis. 

Material for inclusion in these collections should be sent to the Com- 
putation Laboratory, National Bureau of Standards, Washington 25, D. C., 
for the attention of J. H. WEGSTEIN. 


16. Remington Rand, Prefab Tab. A routine is described which generates a 
routine for turning UNIVAC into the logical equivalent of a fictitious 
punch card tabulator. 

17. Remington Rand, Preliminary Exposition, UNIVAC Short Code. 

18. Willow Run Research Center, University of Michigan, Willow Run 
Research Center Memoranda. The following subjects are included : 

MIDAC characteristics 

A list of problems and machines with which they were solved 
Programming on the Whirlwind Computer 

Indoctrination of potential coders 

A suggested method of processing problems on MIDAC 
MIDAC Logic: word composition 

MIDAC magnetic drum control equipment 

Preliminary coding manua! for MIDAC 

19. Analysis Laboratory, Calif. Inst. of Technology, Bibliography of articles 
on computing machines published prior to June, 1949. 

20. Digital Computer Laboratory, Mass. Inst. of Technology, A Short 
Guide to Coding, May, 1951. 

21. Digital Computer Laboratory, Mass. Inst. of Technology, Short Guide 

to Coding and Whirlwind I Operation Code. 
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22. Digital Computer Laboratory, Mass. Inst. of Technology, Programming 
for Whirlwind I. This report tells what Whirlwind I is and what it can 
do. It presents essentials and examples of Whirlwind I programming and 
suggests various programming techniques. It contains an appendix which 
discusses handling numbers in the machine and the operation code. 

23. National Bureau of Standards Report 1612, April 21, 1952. 

SEAC Operating and Programming Notes, I. 


1. 


IA MNP WN? 


Operation of SEAC Control Panel as of October, 1951 


. Non-mathematical Routines 

. Conversion Subroutines for Integers 

. Decimal to Binary and Binary to Decimal Subroutines 
. Subroutine for 2? and e? 

. Subroutines for VN, VN, VN 

. Subroutine for Sin X and Cos X 


24. National Bureau of Standards Report 1807, July 22, 1952. 
a Operating and Programming Notes, II. 


. Subroutine for tan a/b 


. Hexadecimal-decimal Converter Table 
. Four-Hexi Converter Table 
. Subroutine for Decimal to Binary Conversion of a Double Precision 


Number with Fixed Binary Point 


. Subroutine for i-th Differences, i=1, 2, ..., 6 
. Three-address System and Base Order 
. Preparation of Pure Hexadecimal Coding for Punched-Card Com- 


position 


25. National Bureau of Standards Report 1873, August 22, 1952. 
SEAC Operating and Programming Notes, III. 


16. 
17. 


26 


Magnetic Tape to Wire Transfer and Check Routine 
Subroutines for Basic Arithmetic Operations and Square Root in 
Floating Decimal Form 


. General Dump Routine (512 or 1024 Word Memory) 
. Memory Decomposition Routine (512 or 1024 Words) 


. Flow Charts for Magnetic Tape Checking Routines 


National Bureau of Standards Report 1917, September 12, 1952. 


SEAC Operating and Programming Notes, IV. 


21. 
22. 


Composition Routine (SEBBE) 
Basic Arithmetic Operations I, Floating Binary Point, Single Pre- 
cision Numbers 


. Subroutine for Cube Root: Single Precision, Floating Binary Point 
. Subroutine for sin z, cos 2, 3 = x + ty. 


. Corrections and Changes to Technical Memoranda Nos. 1-24. 


27. Rut B. Horcan, SWAC Coding Manual, National Bureau of 
Standards Report 2047, November 4, 1952, 37 p. 


This coding guide gives a good introduction to the art of coding for the 
SWAC, although it is incomplete in certain respects. For example, the rela- 
tive time to obey each command is not given so it is not possible to minimize 
or even estimate the time required to do a calculation with the SWAC. The 
number of binary digits in a word is not stated explicitly, nor whether the 
modulus or complement convention is used. 
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In Chapters 2 and 3 the input-output procedures and possibilities are 
described in detail. Chapter 4 is devoted to a method of verifying the cor- 
rectness of the store and the coding precautions which have to be taken to 
enable this to be done readily. The following two chapters describe how 
commands can be modified and tallies kept. 

The last chapter describes the use of sub-routines and the necessary con- 
ventions although it is not actually stated how a sub-routine is entered from 
a program. 

D. J. WHEELER 
University of Illinois 
Urbana, Illinois 


BIBLIOGRAPHY Z 


1038. ALBERT A. GIRLACH, “Test pulse generator for digital computers,” 
Electronics, v. 25, Nov. 1952, p. 158-160. 


This paper gives a description of a pulse generator which was built using 
multivibrator and gate circuits. It can be triggered either internally or ex- 
ternally. The outputs are a clock output plus two pulse outputs and their 
complements. The width and amplitude of each output can be varied inde- 
pendently. A block diagram and complete circuit diagrams are given. 


ERNEST F. AINSWORTH 
NBSEC 


1039. S. E. Giuck, “The Electronic Discrete Variable Computer,” Elec- 
trical Engineering, v. 72, Feb. 1953, p. 159-162. 


A careful description of the EDVAC is given by the author. This ma- 
chine has been operating successfully since the Spring of 1952 at the Ballistic 
Research Laboratories, Aberdeen Proving Ground, Md. 

The reader familiar with the SEAC will find much similarity between 
the two machines. 

IDA RHODES 
NBS 


1040. INTERNATIONAL BUSINESS MACHINES CORPORATION, Industrial Com- 
putation Seminar, edited by Cuthbert C. Hurd, Dec. 1949, 173 pages. 
22 X 28.5 cm. 


This is one of a series of seminars that was held for research engineers 
and scientists who have active interests in the mathematical solution of 
physical problems by the application of punched card techniques to the 
computing art. The following articles were discussed : 


1. ‘The future of high-speed computing,”’ by JoHN von NEUMANN. 

2. “Some methods of solving hyperbolic and parabolic partial differ- 
ential equations,’’ by RICHARD W. HAMMING. 

3. ‘“‘Numerical solution of partial differential equations,” by EVERETT 
C. YOWELL. 

4. “An eigenvalue problem of the Laplace operator,” by Harry H. 

HUMMEL. 
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10. 
11. 


12. 
13. 


14. 
15. 
16. 
17. 
18. 
19. 
20. 
21. 
23. 
24. 


25. 
26. 


27. 
28. 
29. 
30. 


31. 
32. 
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. “A numerical solution for systems of linear differential equations 


occurring in problems of structures,” by Pau E. Biscu. 


. ‘Matrix methods,” by Kaiser S. Kunz. 
. “Inversion of an alternant matrix,’’ by BoNALYN A. LUCKEY. 
. “Matrix multiplication on the IBM card-programmed electronic 


calculator,”’ by JoHNn P. KELLy. 


. “‘Machine methods for finding characteristic roots of a matrix,” by 
FRANz L. ALT. 
“Solution of simultaneous linear algebraic equations using the IBM 


Type 604 Electronic Calculating Punch,” by Joun Lowe. 
“Rational approximation in high-speed computing,” by CEcIL 
HASTINGS, JR. 

“The construction of tables,” by PAuL HERGET. 

“A description of several optimum interval tables,”’ by Stuart L. 
CROSSMAN. 

“Table interpolation employing the IBM Type 604 Electronic 
Calculating Punch,” by EVERETT KIMBALL, Jr. 

“An algorithm for fitting a polynomial through m given points,” by 
F. N. FRENKIEL and H. PoLacHex. 

“The Monte Carlo Method and its application,’”’ by MARK Kac and 
M. D. DonsKER. 

“‘A punched card application of the Monte Carlo Method,” by P. C. 
JOHNSON and F. C. UFFELMAN. 

““A Monte Carlo Method of solving Laplace’s equation,” by 
EvERETT C. YOWELL. 

“Further remarks on stochastic methods in quantum mechanics,” 
by GILBERT W. KInNc. 

“Standard methods of analyzing data,” by JoHN W. TUKEY. 

“The application of machine methods to analysis of variance and 
multiple regression,”’ by ROBERT J. MONROE. 

“Examples of enumeration statistics,” by W. WAYNE COULTER. 
“Transforming theodolite data,” by HENRY SCHUTZBERGER. 
“Minimum volume calculations with many operations on the IBM 
Type 604 Electronic Calculating Punch,” by WiLL1aM D. BELL. 
“Transition from problem to card program,” by GREGORY J. TOBEN. 
“Best starting values for an iterative process of taking roots,’ by 
PRESTON C. HAMMER. 

“Improvement in the convergence of methods of successive approxi- 
mations,” by L. RICHARD TURNER. 

“Single order reduction of a complex matrix,’’ by RANDALL E. 
PORTER. 

“Simplification of statistical computations as adapted to a punched 
card service bureau,” by W. T. SourHwortu and J. E. BACHELDER. 
“Forms of analysis for either measurement or enumeration data 
amenable to machine methods,”’ by A. E. BRANDT. 


“Remarks on the IBM Relay Calculator,” by M. Lorxin. 


“An improved punched card method for crystal structure factor 
calculations,” by M. D. Gres. 
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33. ‘The calculation of complex hypergeometric functions with the 
IBM Type 602-A Calculating Punch,” by Harvey GELLMAN. 
34. “The calculation of roots of complex polynomials using the IBM 
Type 602-A Calculating Punch,” by Jonn Lowe. 
35. “‘Practical inversion of matrices of high order,” by WiLi1am D. 
GUTSHALL. 


Several of the papers that were discussed dealt with specific problems as 
encountered in actual practice and their particular solution and are useful 
only in that they illustrate the type of approach that the authors used. 
Some of the other articles, such as no. 3 by E. C. Yowell, describe methods 
which are rather general in approach and could be applied to a large number 
of problems involving solution of partial differential equations. The article 
by Donsker and Kac wherein they discuss the Monte Carlo Method and its 
application is also rather general and so could be applied to numerous prob- 
lems. Finally the articles by Kunz on matrix methods and Tukey on the 
standard methods of analyzing data are more theoretical and are excellent 
summaries but contain no mention of any direct application to punch card 
methods. 

H. BREMER 
NBSCL 


1041. HERBERT KopPEL, “Digital computer plays NIM,” Electronics, v. 25, 
Nov. 1952, p. 155-157. 


Engineers of the W. L. Maxon Corporation have constructed an auto- 
matic machine for playing NIM with a human opponent. It can be adjusted 
so that the human can win if he plays a perfect game or so that he can never 
win. Digital computer techniques such as gates, binary counters, etc., were 
used in its construction. The method of operation and a block diagram are 
given but no circuitry details are shown. 

ERNEST F. AINSWORTH 
NBSEC 


1042. W. H. MacWiu1ams, Jr., ‘“Computers—past, present, and future,” 
Electrical Engineering, v. 72, Feb. 1953, p. 116-121. 


The author gives an exceedingly clear, yet non-technical, exposition of 
many computing devices, both analogue and digital, which humanity has 
been utilizing for the past three centuries. This informative article should 
appeal to every member of that rapidly growing circle of engineers and 
mathematicians, who—either as constructors or users—are closely con- 
nected with the art of computation. 

The reviewer is disappointed that the author’s glimpse into the future 
failed to perceive a much needed system of electronic sorter-collators, whose 
lack is felt by agencies handling huge masses of data. 

IpA RHODES 
NBSCL 


1043. OFFICE OF NAVAL RESEARCH, Digital Computer Newsletter, v. 5, 
January 1953, 8 p. 
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The contents are as follows: 


owvaonn 


_ 


1. 
2. 


1. Naval Proving Ground Calculators 
2. Whirlwind I 

3. 
4 
5 


Moore School Automatic Computer 


. The SWAC 
. Aberdeen Proving Ground Computers: 


The ORDVAC 
The EDVAC 
The ENIAC 
The BELL 
IBM-CPC 


. The Circle Computer 

. The Jacobs Instrument Company Computer (JAINCOMP) 
. The ELECOM Computers 

. University of Illinois Computer (ILLIAC) 

. Hughes Aircraft Company Computer 


Data Processing and Conversion Equipment 
Flying Typewriter 
The Charactron 


List of Computing Services 


1044. OrFICE OF NAVAL RESEARCH, Digital Computer Newsletter, v. 5, 


April 1953, 12 p. 


The contents are as follows: 


1. 
. Whirlwind I 
. Computer Research Corporation Computers 





Naval Proving Ground Calculators 


CADAC 102-A 
CRC 105 
CRC 107 


. Moore School Automatic Computer (MSAC) 

. Air Force Missile Test Center Computer (FLAC) 

. The SEAC 

. The [AS Computer 

. The SWAC 

. The MONROBOT 

. The Circle Computer 

. The Jacobs Instrument Company Computers (JAINCOMPS) 
. Consolidated Electronic Digital Computer Model 30-201 
. The ERA 1103 Computer 

. The Rand Corporation Computer 

. Aberdeen Proving Grounds Computers 


Data Processing and Conversion Equipment 


. TELEDUCER 
. SADIC 
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3. Benson-Lehner Incremental Plotter 
4. Coleman Digitizer 
5. Ferro-Resonant Flip-Flop 
6. Logrinc Automatic Graph Followers 
List of Computing Services 
Computer and Numerical Analysis Courses 
1. Massachusetts Institute of Technology 
2. Computer Research Corporation 


1045. L. PAcKER & W. J. Wray, JR., “Germanium photo-diodes read com- 
puter tapes,’ Electronics, v. 25, Nov. 1952, p. 150-151. 


This article describes a system for reading Teletype paper tapes using 
1N77 photo-diodes. Complete circuit details are given for the amplifier and 
shaper. It operates over a range from one to 35,000 pulses per second. 
Mechanical details of the tape transport are not discussed. 


ERNEST F. AINSWORTH 
NBSEC 


1046. UNIVERSITY OF SIDNEY, Proceedings of Conference on Automatic Com- 
puting Machines, held in the Department of Electrical Engineering, 
University of Sidney, Australia, in conjunction with the Commonwealth 
Scientific and Industrial Research Organization, Aug. 1951, 220 pages. 
20 X 24.8 cm. 


The first of two sessions contained ‘‘An introduction to automatic cal- 
culating machines” and a paper entitled “Automatic digital calculating 
machines” both by D. R. HARTREE. Two other papers discussed the C.S.I. 
R.O. (Commonwealth Scientific and Industrial Research Organization) 
Differential Analyser and the C.S.I.R.O. Radiophysics MK. I Automatic 
Computer. The latter is an acoustic-delay-line-memory digital computer 
with magnetic drum and punched card input-output. In the discussion of 
the use of superfluous binary digits for error-detecting which followed, 
Hartree considered that ‘‘the use of error detecting procedures at each 
operation of the machine was a counsel of despair.” He felt that experience 
has shown that computers do not go wrong often enough to warrant this. 

In the second session, Hartree gave an introduction to programming 
using the EDSAC for illustration; also he presented a paper on numerical 
methods used with automatic calculating machines. The latter included a 
summary of various methods of determining roots of polynomial equations 
with automatic machines. T. PEARCEY presented papers on programming 
for the C.S.I.R.O. digital machine and for punched-card machines and on 
the functional design of an automatic computer. Some of the other papers 
presented were: (1) ‘Some analogue computing devices,” (2) “Digital- 
analogue conversions,’ (3) ‘‘An analogue computer to solve polynomial 
equations with real coefficients,” and (4) ‘‘Some new developments in equip- 
ment for high-speed digital machines.’’ The last paper dealt with a high- 
speed magnetic switching device, a single electron tube scale-of-ten numeral- 
displaying counter, a single tube which combines the functional operations 
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of a bi-stable element and gate, giving several applications and details of a 
magnetic-drum digital storage system. 

J. H. WEGsTEIN 
NBSCL 


1047. B. G. WELBY, “‘Intermittent-feed computer tape reader,’’ Electronics, 
v. 26, Feb. 1953, p. 115-117. 


A tape reader is described which will operate at speeds up to 200 char- 
acters per second with teletype tape. Each character consists of as many as 
seven holes punched in the paper tape. Reading is accomplished by a 
photoelectric system, and no reading pins or contacts are required. The tape 
is friction driven by rollers which can be controlled intermittently. The 
maximum reading speed under intermittent control is not stated but ap- 
pears to be about 100 characters per second. Reading causes no visible wear 
to the tape after 10,000 passes. 

5. &. Pea 
NBSCL 


NEws 


ATEE-IRE-ACM.—The joint AIEE-IRE-ACM Computer Conference Committee 
sponsored a Western Computer Conference at Los Angeles, California, on February 4, 5, 
and 6. The program for the meeting was as follows: 

Feb. 4, 1953, 8:30 a.m. 


Registration Ballroom Floor 


10:30—11:30 a.m. 


Opening ceremonies P. L. Morton, Chairman, Univ. of Calif., 
Berkeley 
Keynote addresses : 
The impact of computer development on S. Ramo, Vice-President for Operations, 


the training and utilization of engineers Hughes Aircraft Co., Culver City, Calif. 
Factors influencing the effective use of R. D. Huntoon, Chief, Corona Labora- 
computers tories, NBS Corona, Calif. 
Lunch G. D. McCann, Toastmaster, Calif. Inst. of 
Tech., Pasadena 
The scientific manpower problem L. A. DuBrincz, Pres., Calif. Inst. of Tech., 
Pasadena 
2:00-4:30 p.m. 
Session I, Commercial applications E. C. Netson. Chairman, Hughes Aircraft 


Co., Culver City, Calif. 
Commerical applications—The implica- J. L. McPHERsON, Bureau of the Census 
tions of Census experience 
Payroll accounting with ELECOM 120 R. F. SHaw, Electronic Computer Division, 
Underwood Corp. 
Automatic data-processing in larger man- M. W. SALVESON and R. G. CANNING, Univ. 


ufacturing plants of Calif., Los Angeles 
The data-processing requirements of the E. E. StickeLt, Bureau of Old-Age and 
Social Security problem Survivors Insurance 


The processing of information-containing G. W. Brown and L. N. RIpENovr, Inter- 
documents national Telemeter Corp. 
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Feb. 5, 1953, 9:00-11:30 a.m. 
Session II, Applications to aircraft and C.SrtranG, Chairman, Douglas Aircraft Co., 


missile design 

Landing gear simulation using a differ- 
ential analyzer 

The equivalent circuits of shells used in 
airframe construction 

Analog-digital techniques in autopilot 
design 

Applications of computers to aircraft dy- 
namics problems 


Santa Monica, Calif. 
D. W. Drake, Lockheed Aircraft Corp. 


R. H. MAcNEAL, Calif. Inst. of Tech. 


R. L. JoHNsoN and W. T. Hunter, Douglas 
Aircraft Co. 

D. Dit, B. Hatt, R. RutHraurr, Douglas 
Aircraft Co. 


12:00-1:30 p.m. 


Lunch 


Luncheon Address: 
New equations for management 


R. G. CANNING, Toastmaster, Univ. of Calif., 
Los Angeles 


J. E. Hosson, Director, Stanford Research 


Institute, Palo Alto 


2:00-4:30 p.m. 


Session III, Panel discussion 
An evaluation of analog and digital com- 
puters 


G. D. McCann, Moderator, Calif. Inst. of 
Tech. 


Panel: J. L. BARNEs, Assoc. Director of Electro-Mechanical Engineering Dept., North 
American Aviation, Inc. 
L. RmEnour, Vice-President, International Telemeter Corp. 
F. STEELE, Vice-President in charge of Engineering Digital Control Systems, La 











Jolla, Calif. 


A. W. VANCE, Research Section Head, RCA Laboratory, Princeton, N. J. 


Feb. 6, 1953, 9:00-11:30 a.m. 


Session IV, New developments in digital 
computer equipment 
The snapping dipoles of ferro-electrics as 
a memory element for digital computers 
Magnetic reproducer and printer 


An improved cathode ray tube storage 
system 

Nonlinear resistors in logical switching 
circuits 


H. D. HusKkey, Chairman, Wayne Univ. 


C. F. Putvari, The Catholic Univ. of 
America 

J. C. Sms, Jr., Eckert-Mauchly Div., Rem- 
ington Rand, Inc. 

R. THORENSEN, NBSINA 


F. A. ScHwertz, and R. T. STEINBACK, 
Mellon Institute 


1:30-4:00 p.m. 


Session V, New developments in analog com- 
puting equipment 
New laboratory for a three dimensional 
guided missile analysis 
A new concept in analog computers 


C. H. Wits, Chairman, Calif., Inst. of 
Tech., Pasadena 

L. Baver, Head, Project Cyclone, Reeves 
Instrument Corp. 

L. Cann, Chief Computer Engineer, Special 
Products Dept., Beckman Instruments, 
Inc. 











204 AUTOMATIC COMPUTING MACHINERY 


A magnetically coupled low-cost high- A. J. WINTER, Research Lab., Supervisor, 
speed shaft position digitizer Telecomputing Corp. 

Solution of partial differential equations R.M. Howe, Univ. of Michigan 
by difference methods using the elec- 
tronic differential analyzer 

A syncro-operated differential analyzer A. NorpDsIEcK, Univ. of Illinois 


Following is a brief description of some of the main points brought out in the session on 
commercial applications. Mr. McPherson’s talk outlined conclusions drawn from Census 
experience with the UNIVAC concerning the commercial applicability of the electronic 
digital computers. It is possible for them to exceed in efficiency any other available tool for 
many commercial purposes. They are not commercially ideal because their arithmetic power 
exceeds their input-output power. It is hoped that these machines will evolve into a powerful 
aid to business problems similar to the evolution of punched card equipment. The ELECOM 
120 described by Mr. Shaw, can process the weekly payroll for 4,000 to 5,000 employees at 
the rate of about 30 seconds per employee including typing of checks and statements. It isa 
moderately priced decimal computer. In the discussion by Messrs. Salveson and Canning, 
it is mentioned that the ONR is sponsoring a project on “Production planning and schedul- 
ing’”’ which includes the systems design specifications. Mr. Stickell discusses the possibilities 
and limitations of card-to-tape and tape-to-card equipment for sorting data in electronic 
media and the present lack of equipment for high-speed random access to data stored outside 
of a machine. In the last talk of Session I, Messrs. Brown and Ridenour state that more than 
10,000 documents of all sorts can be processed per hour by a single machine which makes no 
demands on the size, shape, thickness, flatness, or degree of preservation of the documents 
being handled. 

In Session II, “Applications to aircraft and missile design,” Mr. Drake states that the 
computer simulation includes velocity square orifice damping, a nonlinear oleo air spring, 
wheel spin-up, slop in gear joints, and a broad range of airplane weights, effective wing lift, 
and landing speeds. The study suggests that this method of investigating actual landing 
gears is practical. Mr. MacNeal describes the development of analog computer techniques 
that should be heipful in solving practical problems in connection with aircraft fuselage and 
wing design. In the third paper of this session the roles assigned to digital and analog com- 
putation connected with autopilot design at Douglas are discussed and the associated rea- 
sons for such assignment are given. Messrs. Dill, Hall, and Ruthrauff discuss the basic 
methods for the solution of aircraft dynamics problems on both analog and digital com- 
puters; in this connection actual problems successfully solved on the computing equipment 
used by the Douglas Aircraft Company are outlined. 

Mr. Pulvari begins Session IV with a discussion of ‘The snapping dipoles of ferroelec- 
trics as a memory element of digital computers.’’ The following description of his talk is 
quoted from the program: A sensitive pulse method has been developed for obtaining static 
remanent polarization data for ferroelectric materials. This method has been applied to 
study the effect of pulse length and amplitude, and decay of polarization on ferro-electric 
ceramic materials with fairly large crystalline orientation. Attempts have been made to 
develop electrostatically-induced memory devices using ferroelectric substances as a 
medium for storing information, particularly delay line, matrix, and tape type of memory 
devices ; also a simple counting circuit using ferroelectric condensers as a bistable element has 
been designed. The paper by J. C. Sims, Jr., describes a new process for producing printed 
copy by magnetic fields which can be adapted to computer output printing. After the data 
are recorded on a magnetized surface, the latent magnet image is developed with a magnetic 
ink and transferred to paper and fixed. In the paper by R. Thorensen, it is shown that by 
altering somewhat the mode of operation of the cathode ray tube and by changing the as- 
sociated gating circuitry, a system is obtained which operates successfully even under severe 
spillover signal distortion. Messrs. Schwertz and Steinback show how nonlinear resistors 
which are made by applying printed circuit techniques to such materials as standard plastic 
may be used to replace whole arrays of crystal rectifiers in certain logical switching circuits. 
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Some new developments in analog computing equipment were described in Session V. 
A new simulation laboratory constructed at Project Cyclone was planned to be large enough 
to permit the solution of complex three-dimensional guided missile problems. Mr. Bauer 
stated that it may also be used for other types of problems. In the paper “A new concept in 
analog computers” circuits of a coordinated low-cost analog computer are described and 
evaluated in terms of computing ability. 

Symposium on automatic digital computation.—A symposium on Automatic Digital 
Computation was held at the National Physical Laboratory, Teddington, England on the 
25th, 26th, 27th and 28th of March, 1953. 

About 200 delegates attended the Symposium and of these about twenty came from 
other European countries and five from the U.S.A. The United States delegates were 
GERTRUDE BLANCH and R. J. Stutz, NBS; R. Hutsizer, University of Illinois; J. C. Mc- 
Puerson, I.B.M. World Headquarters, New York; and S. Kaurman, Shell Development 
Company. 

The computing machine developed at the National Physical Laboratory, the ACE 
Pilot Model, was demonstrated at various times while the Symposium was in progress. 

The programme for the Symposium was as follows: 


March 25, 1953, 11:30 a. m.—1:00 p.m. 


Session 1 E. T. Goopwin, Chairman, NPL 
Opening remarks E. C. BuLiarp, Director, NPL 
Address D. R. HARTREE, Cambridge Univ. 


2:00 p.m. to 5:00 p.m. 


Session 2, British Machines F. M. CoLesrook, Chairman, NPL 

ACE Pilot Model J. H. Witxrinson, NPL 

EDSAC M. V. WiLkEs, Cambridge Univ. 

Operating and engineering experience J. M.M. PINKERTON, Messrs. J. Lyons and 
gained with LEO Co. 

MADAM F. C. WiLitaMs, Manchester Univ. 

MOSAIC, the Ministry of Supply Auto- A. W. M. Coomps, Post Office Research 
matic Computer Station 

NICHOLAS N. D. Hit, Ettiotr Bros. 

Advance notes on RASCAL E. J. PETHERICK, RAE 

The T. R. E. High-Speed Digital Com- A.M. UttLey, TRE 
puter 


March 26, 1953, 9:30 a.m.—1:00 p.m. 
Session 3, Programming E. T. Goopwin, Chairman, NPL 
Optimum coding G. G. ALway, NPL 


Micro-programming and the choice of J. G. StrInGER, Cambridge Univ. 
order code 


Conversion routines E. N. Mutcu and S. Git, Cambridge Univ. 
Getting programmes right S. Grt_, Cambridge Univ. 
2:00-5:00 p. m. 
Session 4, Design J. H. WiLxinson, Chairman, NPL 


H. 
Special requirements for commercial or T. R. THompson, Messrs. J. Lyons and Co. 
administrative applications 
W. Davies, NPL 
O. 


Input and output D. 

Echelon storage systems D. O. CLaypEen, NPL 

Serial digital adders for a variable radix of R. TOWNSEND, British Tabulating Machine 
notation Co. 


































Session 5A, The Utilization of Computing 
Machines—I 
Mathematics and computing 


Linear algebra on the Pilot ACE 

The numerical solution of ordinary differ- 
ential equations 

The solution of partial differential equa- 
tions by automatic calculating ma- 
chines 


Session 6A, The Utilization of Computing 
Machines—II. Mathematical tables 
Applications of electronic machines in 
pure mathematics 
‘Monte Carlo’ methods for the iteration 
of linear operators 


Session 5B, Circuitry and Hardware 
Gates and trigger circuits 


Mercury delay line storage 

Applications of magnetostriction delay 
lines 

Cathode ray tube storage 

Memory studies at the National Bureau of 
Standards, Washington, D. C., U.S.A. 


Session 6B, Servicing and Maintenance 
Preventive or curative maintenance 
Experience with marginal checking and 

automatic routining of the EDSAC 
Diagnostic programmes 
Component reliability in the computing 
machine at Manchester University 


Session 7, Medium-Size Digital Computing 
Machines 
The Harwell Computer 
The A.P.E.(X)C. A low cost electronic 
calculator 
The Elliott-N.R.D.C. Computer 401—A 
demonstration of computer engineering 
by packaged unit construction 
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March 27, 1953, 9:30 a.m.—1:00 p.m. 


E. T. Goopwin, Chairman, NPL 

A. VAN WIJNGAARDEN, Mathematical 
Centre, Amsterdam 

J. H. Wirkinson, NPL 

L. Fox and H. H. Rospertson, NPL 


N. E. Hoskin, Manchester University 


2:00 p.m.—5:00 p.m. 


L. Fox, Chairman, NPL 
E. T. Goopwin, NPL 
J. C. P. MiLLer, Cambridge University 


J. H. Curtiss, NBS (presented by G. 
Blanch) 


9:30 a.m.—1.00 p.m. 


E. A. NEwMAN, Chairman, NPL 

W. W. CHANDLER, Post Office Research 
Station 

M. Wricat, NPL 

R. C. Rospsins and R. MILLERsHIP, 
Messrs. Elliott Bros. 

T.KiLBurn, Manchester Univ. 

R. J. Stutz, NBS 


2:00 p.m.—5:00 p.m. 


F. M. CoLEesrook, Chairman, NPL 

E. A. NEwman, NPL 

M. V. WILKEs, M. PuistTeEr, JR. and S. A. 
BaRTEN, Cambridge Univ. 

R. L. GriMsDALeE, Ferranti Ltd. 

A. A. Rostnson, Ferranti Ltd. 


March 28, 1953, 9:30 a.m.—11:00 a.m. 


J. R. WomersLey, NPL 


E. H. Cooke-YarsoroucH, AERE 
A. D. Boorn, Birkbeck College 


W. S. Extiotrr, H. G. CARPENTER and 
A. St. Jounston, Elliott Bros. 
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OTHER AIDS TO COMPUTATION 


OTHER AIDS TO COMPUTATION 
COMPARISON OF AMERICAN ELECTRIC DESK CALCULATORS 


The following list of advantages and disadvantages of the three American 
electric desk calculators is intended to be an impartial survey of non-debat- 
able facts concerning the use of the machines for statistical and technical 
problems. 

Anything common to all three machines—on either side of the ledger— 
is excluded. To save space, each item is entered for the odd machine, that is, 
an advantage for two of them is entered as a disadvantage of the third. 

No attempt is made to assess ruggedness, service, or length of life. 

All three machines have features which, though unique, do not clearly 
belong on one side of the ledger or the other. Examples: 


1. Monroe claims “‘velvet touch,” which refers to all the keys having the 
same pressure to operate. 

2. Marchant features higher cycling speed and continuous mesh mech- 
anism in driving the counter wheels. These two almost exactly compensate 
each other. 

3. Friden offers the possibility of adjusting each wheel of the accumulat- 
ing counter individually. 


The comparison is made between the latest models as of the end of 1952, 
and for each the most complete machine, for which the current prices are 
$775 for Monroe and Friden and $815 for Marchant. The Marchant has 
complete carry-over on all models ; the Monroe has it as a standard item only 
on the full size (automatic) machine ; for the Friden it is not standard on any 
model. It should be pointed out that if the work to be done does not call for 
a full size machine, the prices differ sharply. In all cases, it pays to examine 
the range of models of all three before purchasing. No one of the machines 
stands out for all-purpose operation. 

No comparison is made between these machines and other calculators, 
such as the Olivetti and Remington Rand printing calculators, the Curta 
8-ounce hand operated machine, and the Swedish Facit (the latter formerly 
made here by R. C. Allen). 

Both Monroe and Marchant offer octal calculators in their line. Monroe 
and Friden are about to bring forth electronic machines in team with IBM. 
Friden has a square root model available ($1300). 

It is not possible, without exhaustive research, to evaluate the machines 
in terms of hourly production of written answers. It is probably true that 
for a variety of work (typical in technical applications) the difference in 
production speed is slight. 


Friden advantages: 


1. Zero keys can be raised to lock (positively) any single column. 

2. Rounding (reset to 5) is available in six positions of the accumulating 
dial. 
3. Dials can be cleared manually when locked. 
4. Readily continues a halted division. 
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5. Positive lock for entire keyboard. (When locked, dividend entry is 
suppressed.) : 

6. Has a bell. 

7. Automatic transfer from accumulating to counting dial. 


Friden disadvantages : 


1. Dividend tab key always produces long carriage travel. 

2. Divide keys are arranged so that both must be depressed for normal 
(positive) division. 

3. Plus and minus bars are located outside the shift keys. 

4. The machine is not always ready to use as an adding machine. 

5. Unless some dividend tab key is down, dividend will not enter. 


Monroe advantages: 


1. Lightweight, small, and portable. 

2. Smaller models can be built into a full size machine. 

3. Smaller models (non-automatic multiplication) have a single key to 
zero the entire machine. 

4. Drive shaft is available for hand cranking or minor unjamming. 

5. Factor is entered only once for squaring. 

6. Rounding available (one position only). 

7. Has back transfer from the accumulating dial to the multiplier unit. 

8. A unit counter is available (optional) to tally the number of opera- 
tions. 

9. Single digits of a constant multiplier can be corrected or changed. 

10. Zeros in the multiplier are not entered by key depression. 

11. Hasautomaticcarriage return toselected position after multiplication. 


Monroe disadvantages : 


1. Possible runaway on division (lacks divisor lineup feature). 
2. Lacks positive keyboard lock. 
3. Very easy to depress two keys in one bank at once. 


Marchant advantages: 


1. Automatic return of carriage after division and instant dividend/ 
divisor entry. 

2. Visible keyboard entry dials. 

3. Multiplication can be left to right or right to left. 

4. Multiplication takes place as the multiplier digits are entered; the 
product is available just after the last digit of the multiplier is entered. 

5. Repeat addition or subtraction takes place if the plus bar is depressed 
simultaneously with a multiplier key. 

6. Can have just one key per bank down at a time. 

7. Allows a certain amount of touch control in multiplication with vary- 
ing concavities of the multiplier keys. 


Marchant disadvantages : 


1. Lacks auto dividend entry, but see No. 1 advantage. Generally has to 
be set up for the initial division, but subsequent divisions are programmed 
faster. 
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2. Lacks non-entry of the multiplier. 

3. Lacks provision for constant multiplier. 

4. Lacks dial locks. This makes sums of quotients rather difficult and 
prevents all operations which depend on split dial locks, such as decimal 
accumulation. 

5. Lacks positive keyboard lock. 

6. Algebraic sums of products involves changing controls when the sign 
changes. (Friden and Monroe have a negative accumulative multiplication 
key). 

7. Lacks automatic dial clearance at the start of a multiplication. 


FRED GRUENBERGER 
Numerical Analysis Laboratory 
University of Wisconsin 
Madison, Wisconsin 
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The two dimensional wave equation V?¢ + (w/c)*¢ = 0 is solved for its 
least characteristic frequency by means of an electrical network. The cor- 
responding difference equations lead to a square lattice network with bound- 
ary points at zero potential and with interior points connected to ground 
through ‘“‘negative resistors.’’ These ‘‘negative resistors’ involve an adjust- 
table current source and a galvanometer in a Wheatstone bridge circuit. 
They are current generators which are adjusted until the output current 
equals in absolute value the current that would have flowed in a certain 
resistance connected to ground. However this generated current is opposite 
in direction to the last mentioned current. In each solution, all these nega- 
tive resistances are readjusted until all the galvanometers read zero. In order 
to get a non-zero solution a forcing current is introduced at one point in the 
network and a solution for the difference equation for the above differential 
equation obtained. The frequency w is given an increasing sequence of values 
starting with zero and the amount of forcing current is noted. As w approaches 
the first characteristic frequency, the forcing current approaches zero for a 
given amplitude for the solution. Convergence of the iteration process occurs 
past the first characteristic frequency to a point “‘just below the first anti- 
resonant frequency” and no further. The device will also solve problems in 
forced vibrations. Various vibration problems are illustrated and the possi- 
bility of generalizing the boundary conditions is discussed. 

° F. J. M. 
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NOTES 


152.—AsYMPTOTIC FORMULAS FOR THE HERMITE POLYNOMIALS. Asymp- 
totic formulas for the Hermite polynomials H,,(x) are derived by the method 
of Liouville-Steklov from the integral equation 


exp(— 3x") H,(x) = A, cos(N#x — 3nz) 


4 NA f ” # exp(— 32) H(t) sin(NiLx — ¢]at 
0 


where 


H,(x) = (— 1)” exp(x?) (<) exp(— x?), N=2n+1 
and 
_ [|H,(0)| = n!/(3n)! if m is even 


~_— |N-4H,,’(0)| = N-3(m + 1)!/((n + 1)/2)! if m is odd. 


The method is an iterative process where one obtains the k-th approximation 
by placing the (k—1)-th approximation under the integral sign to obtain a 
series in powers of N~! which is quite accurate when 7 is both large and large 
in comparison with x. SzEG6! and others do not give the actual coefficients 
in the asymptotic series, but content themselves with a proof of the existence 
of the expansions and mention of only the first term or so. Each complete 
iteration in the process becomes progressively more tedious, and yields only 
one more term in the expansion. The writer felt that it would be useful to 
have the explicit expressions for the asymptotic formulas considerably be- 
yond the first term. The expansions which are given below are the result of 
9 iterations for both m even and odd. 

These formulas were tested to calculate Hi9(x) and Hoo(x) for x = 1(1)5 
and they permitted a relative error of about 10~? at x = 1 which increased 
steadily with x to about 10~ at x = 5. 

These formulas can also be used to compute the earlier zeros of H,(x) 
by an iterative process. 

From the definitions of A(x) and B(x) below, we have 


0 = A(x) + N-'tan(N'x — $n) B(x). 
We then employ the iteration formula 
Xign = N-? {(4n + m)ax — arctan(N? A (x,)/B(x,))} 


where m is a suitably chosen integer (positive, negative, or zero). For the 
r-th zero of H,(x) we may use as initial approximation 


N-*(r —4)x meven 


¥o = |N-*(r —1)e odd. 


When applied to Hi, and Hy this process gave a relative error of about 10° 
for the first two zeros and about 10- for the last two. 
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The asymptotic formulas in question may be given as follows. 
exp(— 3x") H,(x)/A, = A(x) cos(Nix — 4nr) + B(x) N-'sin(N'x — 4nr) 
where 


A(x) = Ao(x) + N- Ai(x) +...+ N-* Aa(x) + O(n) 
B(x) = Bo(x) + N— By(x) +...+ N-* Ba(x) + O(n). 


The polynomials A ;(x), B;(x) differ slightly in the two cases: 


Case I n even 


Ao(x) =1 

Ai(x) = — x°/72 + x*/4 

Ao(x) = x/31104 — 11x*/1440 + 19x*/96 

A;3(x) = — x'8/33592320 + 17x*/622080 — 18889x'°/3628800 


+ 1091x*/5760 — 19x*/32 
Ag(x) = x*/67722117120 — 23x”/671846400 + 11153x'*/522547200 
— 177127x"*/43545600 + 14601x*/71680 — 631x*/384 


Bo(x) = x*/6 

B(x) = — x°/1296 + x°/15 — x/4 

B.(x) = x'5/933120 — 7x"/12960 + 901x7/20160 — 19x3/48 
B;(x) = — x/1410877440 + x!7/933120 — 2131x'*/5443200 


+ 4421x°/120960 — 241x'/384 + 19x/32 
Ba(x) = x*7/3656994324480 — 13x%/14108774400 
+ 8503x'°/9405849600 — 400187x'5/1306368000 
+ 21500581x"/638668800 — 7159x7/7680 + 631x*/192 


Case II n odd 
Ao(x) = 1 
Ai(x) = — x8/72 + x°/4 
Ao(x) = x!2/31104 — 11x°/1440 + 19x*/96 — 1/4 
A3(x) = — x!8/33592320 + 17x*/622080 — 18889x'°/3628800 


+ 1111x*/5760 — 21x*/32 
Aa(x) = x*/67722117120 — 23x/671846400 + 11153x'*/522547200 
— 59159x!2/14515200 + 132641x*/645120 — 325x*/192 


+ 21/32 
Bo(x) = 23/6 
By(x) = — x9/1296 + x8/15 — x/4 
Bz(x) = x'5/933120 — 7x"/12960 + 901x7/20160 — 7x*/16 
B;(x) = — x%/1410877440 + x!7/933120 — 2131x"*/5443200 


+ 13333x°/362880 — 1237x°/1920 + 21x/32 
B,(x) = x*7/3656994324480 — 13x*/14108774400 
+ 8503x!°/9405849600 — 400537x'*/1306368000 
+ 7195607x"'/212889600 — 152141x7/161280 + 671x*/192. 


H. E. SAuzer 
NBSCL 


- G. SzEG6, Orthogonal Polynomials. Amer. Math. Soc. Collog. Publications, v. 23, 1939, 
p. 212-213. 
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NOTES 


153.—ANALYTICAL APPROXIMATIONS, [See also Note 143]. The following 
twenty-two approximations are all concerned with the so-called offset 
circle probability function defined by the integral 


(13) 


(14) 
(15) 
(16) 
(17) 
(18) 


(19) 
(20) 
(21) 
(22) 


(23) 


(24) 
(25) 
(26) 
(27) 
(28) 
(29) 


(30) 





g(R,x) = f H+ 15 (px) pdp, 


in which Jo(Z) is the usual Bessel function. Among other properties, 
q(R,x) + g(x,R) = 1 + e-**+2) Jy (Rx). 
The q(R,x) function has been finely tabulated to 6D by the joint 
effort of NBSINA and RAND. 
To better than .0028 over — ~ <x CK ~, 

1 — g(R,x) ; 1 
in a —— oo og , 
ro 1 — g(R,0) (1 + .123x? + .010x*)4 
To better than .0014 over - ~ <x K€ ~, 
q(1,x) = 1 — .393(1 + .093x2 + .007x4)-4. 
To better than .00014 over — © Cx K ~, 
q(i,x) = 1 — .3935(1 + .0968x? + .0047x* + .00028x*)—. 
To better than .0035 over —~ <x KC ow, 
q(2,x) = 1 — .865( 1+.038x? + .004x*)-*. 
To better than .001 over — © Cx CK ~, 
q(2,x) = 1 — .865(1 + .0401x? + .00309x* + .000075x*)—. 
To better than .006 over 0 < y < ~, 
lim . q(R,R - y) = ew’ = 1 
ro 1 —q(R,R) (1 + .015y + .076y? + .040y*)* 
To better than .00037 over0 < y < ~, 
q(.5,.5 + y) = 1 — .1045(1 + .129y + .079y? + .056y*)—*. 
To better than .0007 over 0 < y < =~, 

q(1,1 + y) = 1 — .267(1 + .203y + .079y? + .062y%)—. 
To better than .0009 over 0 < y K ~, 
q(2,2 + y) = 1 — .397(1 + 236y + 06652 + .066y*)-*. 
To better than .0011 over 0 < y 
q(4,4+ y) =1— 45(1+ 22iy oa 
To better than 0013 up Ocy< 
5 
1 


. =. 
lim q(R,R+y) = cS aa t= 1 — (1+.209y+ .061y?+ .062y*)* 
To better than .00011 over 0 < x < 1, 

(1.x) = .6066 + .1500x? — .0238x*. 

To better than .0017 over 0 < x < 2, 

q(2,x) = .135 + .566(x/2)? — .096(x/2)*. 

To better than .0008 over 0 < x < 3, 

q(3,x) = .011 + .231(x/3)? + .654(x/3)* — .329(x/3)°. 

To better than .0019 over 0 < x < 3, 

q(3,x) = [.105 + .930(x/3)? — .282(x/3)*P. 

To better than .0011 over 0 < x < 3, 

q(3,x) = [.105 + .954(x/3)? — .349(x/3)4 + .043(x/3)*P. 

To better than .002 over 0 < x < 4, 

q(4,x) = [.018 + .581(x/4)? + .515(«/4)* — .372(x/4)*. 

To better than .0035 over 0 < y < 3, 

q(3,3 — y) = .568(1 + .157y + .107y? + .017y*)—. 
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NOTES 


(31) To better than .0013 over 0 < y < 4, 
q(4,4 — y) = .551(1 + .187y + .055y? + .051y*)—. 
(32) To better than .0013 over0 C y < ~, 


~~ 2 B 
: ne = —— g—t J; = —_— 
lim q(R,R—-y)= J =e dt * 5005-4 .0619*-+.062y)" 
(33) To better than .0004 over 0 < R < 1, 
q(R,R) = 1 — .4921R? + .3212R* — .0966R°*. 
(34) To better than .0011 over 1 ¢ R < ~, 
q(R,R) = (R + .183)/(2R — .388). 





Cecit HAsTINGs, JR. 
James P. Wonca, Jr. 
The RAND Corporation 
1700 Main Street 
Santa Monica, California 


CORRIGENDA 


V. 7, p. 87, 1. 11 for 6202089 read 62020897 
for 1858477404602617 read 18584774046020617 
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